求解低秩矩陣填充的改進的交替最速下降法

胡劍峰

(海南師范大學 數學與統計學院,海南 海口 571158)

0 引言

隨著科技的高速發展，人類社會已逐步從信息時代進入大數據時代。對大數據的采集、分析與處理在當今的社會生活與科學研究中占據著越來越重要的地位。龐大的數據量往往會為數據采集、存儲、傳輸和處理帶來巨大的困難。幸運地是眾多實際問題中的數據之間往往存在著某種關聯性，一般構成一個低秩或近似低秩矩陣。利用矩陣的低秩性，通過對低維觀測量的某種非線性運算來重構原始矩陣，稱為矩陣恢復，它是壓縮感知的推廣。

矩陣填充是矩陣恢復的一種特殊情形，指的是若只能采樣到原低秩矩陣的部分元素，而其它元素無法得到或在存儲和傳輸過程中丟失，需要由已知的部分元素將其余空缺的元素合理準確地填充恢復出整個原始矩陣。矩陣填充在推薦系統、信號處理、系統識別、醫學成像及機器學習等領域有著廣泛應用[1]。

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

(1)

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

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

(2)

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

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

為克服該計算瓶緊，一些學者考慮如下流形上的優化

(3)

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

另外一類算法則避開SVD的計算，利用矩陣分解Z=XY，其中X∈Rm×r，Y∈Rr×n求解如下非凸優化

(4)

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

1 交替最速下降法

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

(5)

由于對(5)中的最小二乘子問題的精確求解較為復雜，因此文[18]采用精確線搜索的最速下降法對其進行一次迭代近似求解，提出了交替的最速下降法(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)

步長txk和tyk按如下精確線搜索計算得到

(8)

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

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

2 改進的交替最速下降法

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

將(8)中的最小化問題改寫為

(9)

顯然，若將(9)中的變量由數量矩陣形式tI替換為對角矩陣形式，則所得最優值更小，從而使得迭代下降更多。 因此，本文提出如下改進的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。對上述一元凸二次優化可直接求解。

<(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)

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

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

定理1由算法2.1產生的序列{(Xk,Yk)}的每個極限點都是的(4)穩定點。

定理2假設算法2.2產生的每個迭代對(Xk,Yk)都是滿秩的，若其極限點也是滿秩的，則該極限點為(4)的穩定點。

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

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

(15)

(16)

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

(17)

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

(18)

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

(19)

因此

=0

即證。

算法2.1 改進的交替最速下降法(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 改進的調比的交替最速下降法(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 數值結果

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

3.1 隨機矩陣填充

本小節對構造的隨機測試問題進行數值試驗。隨機生成矩陣Z0=XY，其中X∈Rm×r，Y∈Rr×n且元素獨立服從標準正態分布。按均勻分布在Z0中隨機采樣p個元素，對應下標集為Ω。dr=r(m+n-r)是秩為r的m×n矩陣的自由度。由于矩陣的低相干性較難驗證，因而由其確定采樣數p就較困難。然而實際經驗表明[12,14,18]， 一般只需p/dr比1適當地大一些，就能很好地恢復原矩陣。測試問題分為如下三組

* 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列出了各算法求解的相對誤差(reler)、迭代數(itns)和計算時間(time)，其中時間單位為秒，用tic，toc計時，reler=‖Zout-Z0‖F/‖Z0‖F。所有結果均為計算十個隨機問題的平均值。

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

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

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

表1～表3表明對每個隨機試驗問題，所測試算法的相對誤差都可達到10-5的精度，且IASD和IScaledASD的迭代數和計算時間均分別少于ASD和ScaledASD，而IScaledASD又是其中最少的。此外還可發現，當r≤100,p/dr≤3.5時，IASD也超過了ScaledASD。進一步，表4說明了IASD，IScaledASD的總迭代數分別比ASD，ScaledASD減少了15.2%和14.5%，而總計算時間則分別減少了11.3%和10.7%。

表4 總迭代數和計算時間的統計及比較

3.2 圖像修復

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

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

表5 Boat, Barbara和Lena圖像修復的數值結果

4 結論

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