求解低秩矩陣填充的改進(jìn)的交替最速下降法

胡劍峰

(海南師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,海南 ?？?571158)

0 引言

隨著科技的高速發(fā)展，人類(lèi)社會(huì)已逐步從信息時(shí)代進(jìn)入大數(shù)據(jù)時(shí)代。對(duì)大數(shù)據(jù)的采集、分析與處理在當(dāng)今的社會(huì)生活與科學(xué)研究中占據(jù)著越來(lái)越重要的地位。龐大的數(shù)據(jù)量往往會(huì)為數(shù)據(jù)采集、存儲(chǔ)、傳輸和處理帶來(lái)巨大的困難。幸運(yùn)地是眾多實(shí)際問(wèn)題中的數(shù)據(jù)之間往往存在著某種關(guān)聯(lián)性，一般構(gòu)成一個(gè)低秩或近似低秩矩陣。利用矩陣的低秩性，通過(guò)對(duì)低維觀測(cè)量的某種非線性運(yùn)算來(lái)重構(gòu)原始矩陣，稱(chēng)為矩陣恢復(fù)，它是壓縮感知的推廣。

矩陣填充是矩陣恢復(fù)的一種特殊情形，指的是若只能采樣到原低秩矩陣的部分元素，而其它元素?zé)o法得到或在存儲(chǔ)和傳輸過(guò)程中丟失，需要由已知的部分元素將其余空缺的元素合理準(zhǔn)確地填充恢復(fù)出整個(gè)原始矩陣。矩陣填充在推薦系統(tǒng)、信號(hào)處理、系統(tǒng)識(shí)別、醫(yī)學(xué)成像及機(jī)器學(xué)習(xí)等領(lǐng)域有著廣泛應(yīng)用[1]。

設(shè)原始矩陣Z0∈Rm×n，rank(Z0)=r?m,n,Ω為采樣到的部分元素的下標(biāo)集，則可期望通過(guò)求解如下秩最小化問(wèn)題恢復(fù)原矩陣

(1)

然而(1)是NP-Hard的[2,3]。由于rank(Z)等于Z的非零奇異值的個(gè)數(shù)，而核范數(shù)‖Z‖*為Z的奇異值之和，且‖Z‖*為rank(Z)在單位球{Z:‖Z‖*≤1}上的凸包絡(luò)，因此將(1)松弛為核范數(shù)最小化問(wèn)題

min ‖Z‖*
s.t.PΩ(Z)=PΩ(Z0)

(2)

其中PΩ(Z)表示Ω上的正交投影算子，即若(i,j)∈Ω，PΩ(Z)的第(i,j)位置上的元素為Zij，否則為0。Candès和Recht證明了當(dāng)原矩陣具有低相干性且采樣數(shù)滿足一定條件時(shí)，(1)和(2)是等價(jià)的，即有相同的唯一最優(yōu)解且該最優(yōu)解以極大的概率等于原矩陣Z0[4]。其后，Candès和陶哲軒又進(jìn)一步改進(jìn)了其條件[5]。

(2)是一個(gè)凸優(yōu)化問(wèn)題，可求得全局最優(yōu)解。然而對(duì)于核范數(shù)最小問(wèn)題的求解，通常每次迭代都需要計(jì)算一個(gè)m×n矩陣的(部分)奇異值分解(SVD)，如奇異值閾值算法(SVT[6,7])、加速鄰近梯度法(APG[8])、非精確增廣拉格朗日乘子法(IALM[9])、不動(dòng)點(diǎn)算法(FPCA[10])、交替方向法(ADM[11])以及迭代硬閾值算法(IHT[12,13])等。一般m,n是可比的,即m=αn,其中α>0為比例常數(shù)，則SVD的計(jì)算量為O(n3)。當(dāng)m,n很大時(shí)，其計(jì)算量太大，因而制約了上述算法在實(shí)際中的應(yīng)用，尤其對(duì)于需要實(shí)時(shí)求解的問(wèn)題。

為克服該計(jì)算瓶緊，一些學(xué)者考慮如下流形上的優(yōu)化

(3)

其中‖·‖F(xiàn)表示矩陣的Frobenius范數(shù)，Mr是嵌入在Rm×n的(m+n-r)r維光滑子流形。文[14]提出了基于黎曼流形切梯度的非線性共軛梯度法(LRGeomCG)。文[15]提出了基于格拉斯曼流形子空間校正的調(diào)比的梯度法(ScGrassMC)。由于利用了流形結(jié)構(gòu)特點(diǎn)，LRGeomCG和ScGrassMC每次迭代分別計(jì)算一個(gè)2r階和r階矩陣的SVD，計(jì)算量為O(r3)。當(dāng)r?m,n時(shí)，極大地提高了計(jì)算效率。LRGeomCG和ScGrassMC可得到較高精度的解。

另外一類(lèi)算法則避開(kāi)SVD的計(jì)算，利用矩陣分解Z=XY，其中X∈Rm×r，Y∈Rr×n求解如下非凸優(yōu)化

(4)

交替最小化策略簡(jiǎn)單有效, 是求解大規(guī)模矩陣分解問(wèn)題最流行的方法之一[1]。采用交替最小化策略求解(4)的算法主要有：基于PowerFactorization的交替最小化算法(PF[16])，低秩矩陣擬合算法(LMaFit[17]) 和交替最速下降算法及其調(diào)比的變種(ASD和ScaledASD[18])。盡管這些算法目前理論上只能確保收斂到(4)的穩(wěn)定點(diǎn)，但其數(shù)值計(jì)算很有效，尤以最近提出的ScaledASD算法表現(xiàn)更優(yōu)。特別對(duì)大規(guī)模問(wèn)題求解且精度要求中等(如10-5)時(shí)，ScaledASD的計(jì)算效率超過(guò)了LRGeomCG和ScGrassMC[18]。本文在ASD和ScaledASD算法的基礎(chǔ)上, 采用分離地精確線搜索代替原有的精確線搜索，使得每次迭代的目標(biāo)函數(shù)值下降更多而不改變每次迭代的計(jì)算量，從而可進(jìn)一步提高計(jì)算效率。

1 交替最速下降法

求解(4)的PF算法的迭代為

(5)

由于對(duì)(5)中的最小二乘子問(wèn)題的精確求解較為復(fù)雜，因此文[18]采用精確線搜索的最速下降法對(duì)其進(jìn)行一次迭代近似求解，提出了交替的最速下降法(ASD)，迭代如下

Xk+1=Xk-txkfYk(Xk)
Yk+1=Yk-tykfXk+1(Yk)

(6)

fY(X)=-(PΩ(Z0)-PΩ(XY))YT
fX(Y)=-XT(PΩ(Z0)-PΩ(XY))

(7)

步長(zhǎng)txk和tyk按如下精確線搜索計(jì)算得到

(8)

易知(8)中的函數(shù)均為一元凸二次函數(shù)，其最優(yōu)解可直接得到，即

ASD算法每次迭代簡(jiǎn)單且計(jì)算量很小，只需8|Ω|r次浮點(diǎn)運(yùn)算，其中|Ω|表示集合Ω中元素的個(gè)數(shù)。為了加速ASD算法的收斂，提高計(jì)算效率，文[18]又接著提出了調(diào)比的交替最速下降法(ScaledASD)，即將(6)中的負(fù)梯度方向-fYk(Xk)和-fXk+1(Yk)分別替換為類(lèi)似牛頓方向的調(diào)比的負(fù)梯度方向和并同樣采用精確線搜索計(jì)算步長(zhǎng)。ScaledASD算法的單步迭代計(jì)算量只比ASD算法多4(m+n)γ次浮點(diǎn)運(yùn)算。

2 改進(jìn)的交替最速下降法

為了敘述方便，先對(duì)一些記號(hào)加以說(shuō)明。I表示單位矩陣，其維數(shù)符合上下文。diag(x)表示以向量x為對(duì)角元的對(duì)角矩陣。<.,.>表示向量或矩陣的內(nèi)積。xi表示向量x的第i個(gè)分量。Ai.和A.j分別表示A的第i行和第j列,‖·‖表示向量的2-范數(shù)。

將(8)中的最小化問(wèn)題改寫(xiě)為

(9)

顯然，若將(9)中的變量由數(shù)量矩陣形式tI替換為對(duì)角矩陣形式，則所得最優(yōu)值更小，從而使得迭代下降更多。 因此，本文提出如下改進(jìn)的ASD迭代(IASD)

Xk+1=Xk-diag(αk)fYk(Xk)
Yk+1=Yk-fXk+1(Yk)diag(βk)

(10)

其中

(11)

t(PΩ(fYk(Xk)Yk))i.‖2

t(PΩ(Xk+1fXk+1(Yk))).j‖2

(12)

i=1,…,m,j=1,…,n。對(duì)上述一元凸二次優(yōu)化可直接求解。

<(PΩ(Z0-XkYk))i.,(PΩ(fYk(Xk)Yk))i.>

=<(PΩ(Z0-XkYk))i.,(fYk(Xk)Yk)i.>

=<(PΩ(Z0-XkYk))i.,(fYk(Xk))i.Yk>

(13)

(14)

此外，若將上述中的梯度方向替換為其調(diào)比方向且采用相同的分離精確線搜索，則可得到改進(jìn)的ScaledASD算法(IScaledASD)。算法2.1和算法2.2分別描述了IASD和IScaledASD的具體步驟。其中，算法2.1的步驟4中關(guān)于(Xk+1,Yk)處的殘差可按如下遞推計(jì)算PΩ(Z0)-PΩ(Xk+1Yk)=PΩ(Z0)-PΩ(XkYk)+diag(αk)PΩ(fYk(Xk)Yk),算法2.2同理。因此，IASD和IScaledASD的單步迭代計(jì)算量分別與ASD和ScaledASD相同。

下面分析算法2.1和算法2.2的收斂性，得到了與ASD和ScaledASD相同的收斂結(jié)論。由于定理2.1的證明可看作是定理2.2的一種特殊情形，因而本文只給出定理2.2的證明過(guò)程。

定理1由算法2.1產(chǎn)生的序列{(Xk,Yk)}的每個(gè)極限點(diǎn)都是的(4)穩(wěn)定點(diǎn)。

定理2假設(shè)算法2.2產(chǎn)生的每個(gè)迭代對(duì)(Xk,Yk)都是滿秩的，若其極限點(diǎn)也是滿秩的，則該極限點(diǎn)為(4)的穩(wěn)定點(diǎn)。

證明設(shè)(X*,Y*)為子列{(Xkq,Ykq)}的極限。由于fX(Y)和fY(X)為連續(xù)函數(shù)，因而只要證明及即可。

由式(13)和(13)同理可得

(15)

(16)

將式(15)和(16)從k=1累加到k=l可得

(17)

由于f(Xl,Yl)≥0，因此有

(18)

式(18)與文[18]中引理4.2的式(10)完全一樣，條件也一樣，因而根據(jù)其接下來(lái)的證明可得

(19)

因此

=0

即證。

算法2.1 改進(jìn)的交替最速下降法(IASD)Input:PΩ(Z0),X0∈Rm×r,Y0∈Rr×nRepeat1.?fYk(Xk)=-(PΩ(Z0)-PΩ(XkYk))YTk2.αki =‖(?fYk(Xk))i.‖2‖(PΩ(?fYk(Xk)Yk))i.‖2,i=1,…,m3.Xk+1=Xk-diag(αk)?fYk(Xk)4.?fXk+1(Yk)=-XTk+1(PΩ(Z0)-PΩ(Xk+1Yk))5.βkj=‖(?fXk+1(Yk)).j‖2‖(PΩ(Xk+1?fXk+1(Yk))).j‖2,j=1,…,n6.Yk+1 =Yk-?fXk+1(Yk)diag(βk)7. k=k+1Until終止條件滿足Output:Z=XkYk

算法2.2 改進(jìn)的調(diào)比的交替最速下降法(IScaledASD)Input:PΩ(Z0),X0∈Rm×r,Y0∈Rr×mRepeat1. ?fYk(Xk)=-(PΩ(Z0)-PΩ(XkYk)))YTk,dxk=-?fYk(Xk)(YkYTk)-12.αki =-<(?fYk(Xk))i.,(dxk)i.>‖(PΩ(dxkYk))i.‖2,i=1,…,m3.Xk+1=Xk+diag(αk)dxk4.?fXk+1(Yk)=-XTk+1(PΩ(Z0)-PΩ(Xk+1Yk)),dyk=-(XTk+1Xk+1)-1?fXk+1(Yk)5.βkj=-<(?fXk+1(Yk)).j,(dyk).j>‖(PΩ(Xk+1dyk).j‖2,j=1,…,n6.Yk+1=Yk+dykdiag(βk)7.k=k+1Until終止條件滿足Output:Z=XkYk

3 數(shù)值結(jié)果

本節(jié)通過(guò)數(shù)值試驗(yàn)對(duì)交替最速下降法及其調(diào)比變種(ASD和ScaledASD)與本文提出的改進(jìn)算法(IASD和IScaledASD)進(jìn)行比較。其中，ASD和ScaledASD的程序是在原文作者個(gè)人主頁(yè)上下載的matlab程序以及兩個(gè)利用稀疏結(jié)構(gòu)的C語(yǔ)言子程序。IASD和IScaledASD的程序是在其基礎(chǔ)上修改的，且同樣利用了稀疏結(jié)構(gòu)采用C語(yǔ)言編寫(xiě)計(jì)算αk,βk及其對(duì)應(yīng)的對(duì)角矩陣乘積的子程序并用mex編譯。所有數(shù)值結(jié)果均是在Windows 7系統(tǒng)，Intel Core i5- 4460 3.2 Ghz處理器、8GB內(nèi)存的個(gè)人電腦的Matlab 2013b上運(yùn)行所得。

3.1 隨機(jī)矩陣填充

本小節(jié)對(duì)構(gòu)造的隨機(jī)測(cè)試問(wèn)題進(jìn)行數(shù)值試驗(yàn)。隨機(jī)生成矩陣Z0=XY，其中X∈Rm×r，Y∈Rr×n且元素獨(dú)立服從標(biāo)準(zhǔn)正態(tài)分布。按均勻分布在Z0中隨機(jī)采樣p個(gè)元素，對(duì)應(yīng)下標(biāo)集為Ω。dr=r(m+n-r)是秩為r的m×n矩陣的自由度。由于矩陣的低相干性較難驗(yàn)證，因而由其確定采樣數(shù)p就較困難。然而實(shí)際經(jīng)驗(yàn)表明[12,14,18]， 一般只需p/dr比1適當(dāng)?shù)卮笠恍?，就能很好地恢?fù)原矩陣。測(cè)試問(wèn)題分為如下三組

* Test set 1:m=n=8000,p/dr=3,r∈{40,60,80,100,120,140,160}；

* Test set 2:r=80,p/dr=3,m=n∈{4000,6000,8000,10000,12000,16000,20000}；

* Test set 3:m=n=8000,r=80,p/dr∈{2,2.5,3,3.5,4,4.5,5}。

表1～表3列出了各算法求解的相對(duì)誤差(reler)、迭代數(shù)(itns)和計(jì)算時(shí)間(time)，其中時(shí)間單位為秒，用tic，toc計(jì)時(shí)，reler=‖Zout-Z0‖F(xiàn)/‖Z0‖F(xiàn)。所有結(jié)果均為計(jì)算十個(gè)隨機(jī)問(wèn)題的平均值。

表1 Test set 1的數(shù)值結(jié)果(m=n=8000,p/dr=3)

表2 Test set 2的數(shù)值結(jié)果(r=80,p/dr=3)

表3 Test set 3的數(shù)值結(jié)果(m=n=8000,r=80)

表1～表3表明對(duì)每個(gè)隨機(jī)試驗(yàn)問(wèn)題，所測(cè)試算法的相對(duì)誤差都可達(dá)到10-5的精度，且IASD和IScaledASD的迭代數(shù)和計(jì)算時(shí)間均分別少于ASD和ScaledASD，而IScaledASD又是其中最少的。此外還可發(fā)現(xiàn)，當(dāng)r≤100,p/dr≤3.5時(shí)，IASD也超過(guò)了ScaledASD。進(jìn)一步，表4說(shuō)明了IASD，IScaledASD的總迭代數(shù)分別比ASD，ScaledASD減少了15.2%和14.5%，而總計(jì)算時(shí)間則分別減少了11.3%和10.7%。

表4 總迭代數(shù)和計(jì)算時(shí)間的統(tǒng)計(jì)及比較

3.2 圖像修復(fù)

本小節(jié)對(duì)圖像修復(fù)問(wèn)題進(jìn)行數(shù)值試驗(yàn)。和文[18]一樣，測(cè)試圖像為三個(gè)標(biāo)準(zhǔn)的512×512像素的灰度圖(Boat, Barbara和Lena)的秩為50的最佳Frobenius范數(shù)逼近的低秩圖像?？紤]如下圖像修復(fù)：按均勻分布隨機(jī)采樣35%的像素，填充和修復(fù)其余65%的像素。此外，終止條件與初始點(diǎn)選取也與文[18]一樣。表5列出了各算法的數(shù)值結(jié)果，其中‘-’表示達(dá)到迭代次數(shù)上限(5000次)時(shí)仍未滿足相對(duì)殘差終止條件而終止算法。此外，以Boat圖像為例，在圖1中顯示了各算法迭代100次時(shí)對(duì)Boat圖像的恢復(fù)效果。由表5和圖1可知，IScaledASD和ScaledASD遠(yuǎn)比IASD和ASD有效，而IScaledASD和IASD分別優(yōu)于 ScaledASD和ASD。

圖1 Boat圖像的重構(gòu)圖(itns=100)

表5 Boat, Barbara和Lena圖像修復(fù)的數(shù)值結(jié)果

4 結(jié)論

本文采用分離地精確線搜索，給出了求解矩陣填充的一種改進(jìn)的交替最速下降法及其調(diào)比變種，即IASD和IScaledASD。與ASD，ScaledASD一樣，IASD和IScaledASD也屬于一階算法，因而當(dāng)問(wèn)題病態(tài)或需要求得高精度解時(shí)，其同樣會(huì)面臨漸進(jìn)收斂變得緩慢的問(wèn)題。對(duì)于此種情形，可先采用IASD或IScaledASD快速求得一個(gè)中等精度的解作為初始解，再采用高階算法，如LRGeomCG和ScGrassMC,接著求解。此外, 本文討論的是秩已知的問(wèn)題。若矩陣的秩事先不知道，可同樣按照文[18]中的建議，基于采樣矩陣的奇異值分解給出原矩陣秩的一個(gè)估計(jì)，然后采用IASD求解。