稀疏重構(gòu)的一種新的加速動量梯度投影法

武高玉，李慧云

（1.河北工業(yè)大學(xué) 理學(xué)院，天津 300401；2.河北工業(yè)大學(xué) 控制科學(xué)與工程學(xué)院，天津 300130）

稀疏重構(gòu)的一種新的加速動量梯度投影法

武高玉1，李慧云2

（1.河北工業(yè)大學(xué) 理學(xué)院，天津 300401；2.河北工業(yè)大學(xué) 控制科學(xué)與工程學(xué)院，天津 300130）

提出了一種新的應(yīng)用于稀疏信號重構(gòu)的加速動量梯度投影法.該方法是把負梯度方向與動量項的凸組合作為搜索方向，步長采取滯后最速下降法（LSD）的步長選取規(guī)則.與加速動量梯度投影法取固定的學(xué)習(xí)速率和動量參數(shù)不同，該方法是動態(tài)地選取動量參數(shù)和步長，從而加速了算法的收斂.數(shù)值試驗表明，與已有求解大規(guī)模l1正則化最小二乘問題的一些方法相比較，本文提出的算法無論是在時間上還是在信號重構(gòu)的質(zhì)量上都是有競爭力的.

動量梯度投影法；滯后最速下降法；l1正則化；信號重構(gòu)；圖像復(fù)原

0 引言

壓縮感知（Compressed Sensing，CS）理論是最近幾年興起的一種新穎的信號采樣恢復(fù)理論.一經(jīng)提出，便在信號處理[1-3]，機器學(xué)習(xí)與模式識別，圖像處理[4-9]等領(lǐng)域受到高度關(guān)注.壓縮感知理論的主要內(nèi)容包括信號的稀疏性、觀測矩陣的設(shè)計和信號的重構(gòu)算法3部分，其中稀疏信號重構(gòu)算法的研究是壓縮感知理論的核心內(nèi)容.重構(gòu)算法的設(shè)計直接關(guān)系到信號的重構(gòu)質(zhì)量，是衡量算法整體性能的重要指標(biāo).

設(shè)x∈Rn是一稀疏信號，將其投影到測量矩陣Ak×n（k

有噪聲的情況下

其中ρ∈Rk表示噪聲.

稀疏重構(gòu)就是將原信號x從測量信號b中完整地重構(gòu)出來，從數(shù)學(xué)意義上講，即為計算欠定線性方程組（1）（或者（2）） 的稀疏解x.一般將其轉(zhuǎn)化為最優(yōu)化問題求解.

本文研究的是受到廣泛關(guān)注的稀疏重構(gòu)的最優(yōu)化模型-l1正則化最小二乘問題，即

在實際應(yīng)用中，問題（3） 的2個難點是非光滑性和大規(guī)模性.最近幾年，人們考慮將不光滑的部分光滑化，為此提出了很多模型和算法.Kim[2]等考慮問題（3） 作為一個凸二次規(guī)劃問題，提出了一種內(nèi)點方法（l1-ls），該方法通過共軛梯度法計算搜索方向.Zhou[4]等將不光滑部分用光滑函數(shù)逼近，然后采取光滑的信賴域方法求解.Figueiredo[1]等將問題（3） 轉(zhuǎn)化為一個界約束凸二次規(guī)劃問題，并通過不同步長的投影梯度法（GPSR_Basic，GPSR_BB）求解；Ma[6]等在經(jīng)典的負梯度方向上添加一個動量項，得到了一種簡單有效的加速動量梯度投影法（MGPSR）.

加速動量梯度投影法是在動量梯度法[10]的基礎(chǔ)上采取投影技術(shù)而得到的.該方法的優(yōu)點是：改善了一階梯度法的收斂性，避免了一階梯度法（如非單調(diào)譜投影梯度法[11-12]）振蕩的缺點，降低了一階梯度法的對局部最優(yōu)解的敏感度；對于大規(guī)模問題，該方法避免二階方法（如牛頓法[13]，擬牛頓法[14]等）計算二階Hessian矩陣所帶來的計算復(fù)雜性，減少了計算時間.但是，加速動量梯度投影法[6]美中不足的是：對學(xué)習(xí)速率和動量參數(shù)比較敏感，通常由實驗者手動挑選，影響了算法的收斂速度.

基于上面的分析，本文提出了一種新的應(yīng)用于稀疏重構(gòu)的加速動量梯度投影法.該方法的搜索方向取為負梯度方向與動量項的凸組合，并采用動態(tài)的動量參數(shù)；在步長的選取上，采取滯后最速下降法[15]的步長選取原則，加速了算法的收斂.數(shù)值試驗表明，與已有求解大規(guī)模l1正則化最小二乘問題的一些方法相比較，本文提出的方法具有更少的計算量和計算時間，且稀疏重構(gòu)的質(zhì)量較好.

本文組織如下：第2節(jié)介紹了一種新的加速動量梯度投影法；第3節(jié)數(shù)值試驗以及試驗結(jié)果的分析；最后是結(jié)束語及參考文獻.

1 算法

本文考慮問題（3） 的等價形式—文獻[1]中的界約束二次規(guī)劃問題（BCQP）：

其中：u=max（x，0）；v=max（-x，0）.故x=u-v，‖x‖1=lTnu+lTnv，ln=[1，1，…，1]T是n維列向量.

問題（4） 寫為標(biāo)準(zhǔn)的BCQP形式：

Figueiredo[1]等通過選取不同的步長的投影梯度法（GPSR_Basic，GPSR_BB）求解式（5），GPSR方法的主要框架為

然而，GPSR選擇的搜索方向是負梯度方向，常常會導(dǎo)致收斂慢和振蕩的現(xiàn)象.

Ma[6]等提出了求解問題（4）的一種簡單有效的加速動量梯度投影法，該方法的搜索方向是通過在負梯度方向上添加一個動量項得到，這樣的搜索方向改善了投影梯度法的收斂速度且減輕了振蕩的現(xiàn)象.加速動量梯度法投影的主要迭代過程為

在這里，η＞0表示學(xué)習(xí)速率，ω∈[0，1]表示動量參數(shù).λk通過精確線搜索得到.

但是，加速動量梯度投影法的學(xué)習(xí)速率和動量參數(shù)是固定的常數(shù).本文提出了一種新的加速動量梯度投影法，該方法實現(xiàn)了動量參數(shù)的自動選取，從而改善了加速動量梯度投影法收斂速度.本文算法的搜索方向是負梯度方向與動量項的凸組合，即

其中ωk∈[0，1]，滿足

其中β＞0是常數(shù).

在確定下降方向之后，選取LSD的步長選取規(guī)則，即

從zk到zk＋1通過投影技術(shù)實現(xiàn)下一次迭代的更新

通過精確線搜索得到λk∈[0，1]，下一次迭代更新為

接下來總結(jié)一下本文算法.

算法1

Step0（初始化）給定z0，Δz0，ωk∈（0，1），令ηk=η0，設(shè)置：k=0；

Step1（計算步）

Step2（線搜索和更新迭代）

計算γk=（Δzk）TBΔzk；若γk=0，那么λk=1；否則，

Step3（更新ηk）

如果γk=0，令ηk＋1=ηmax，否則，由式（14） 計算ηk＋1，

如果滿足終止條件，停止；否則，令ηk=ηk＋1，轉(zhuǎn)Step1.

注：終止條件有如下選擇：

終止條件1：

F*指函數(shù)的最小值或者近似最小值，其中tolP表示允許誤差.

終止條件2：

其中tolP表示允許誤差.

2 數(shù)值試驗

這一節(jié)，用數(shù)值試驗來說明算法的有效性，本節(jié)所有數(shù)值試驗的運行環(huán)境均為MATLAB-R2010a，聯(lián)想（Intel（R）Core（TM）i3-2130 CPU 3.39 GHZ，3.33 GB）.設(shè)置參數(shù)η0=1，tolP=10-3，ηmin=10-30，ηmax=1030，ωk的選取規(guī)則為

表1 幾種算法的CPU時間Tab.1 CPU time of several algorithms

2.1 稀疏信號重構(gòu)試驗

數(shù)值試驗考慮在一個典型的CS環(huán)境中，目標(biāo)是從k觀測信號上，重建一個長度為n的稀疏信號，k< n.在這個例子中，n=4 096，k=1 024，信號x包含160個隨機放置的±1，矩陣A∈Rk×n服從標(biāo)準(zhǔn)的高斯分布，且每行進行標(biāo)準(zhǔn)正交化，b=Ax＋ρ，ρ∈Rk表示方差為σ2的高斯白噪聲，σ2=10-4.參數(shù)τ的選取為：

τ=0.1‖ATb‖∞

把算法1與l1-ls[2]，IST[7]，GPSR_Basic，GPSR_BB，MGPSR[6]，進行比較.為此，運行l(wèi)1-ls算法，并且每一個算法直到達到l1-ls算法得到的目標(biāo)函數(shù)值才終止.

表1是幾類算法應(yīng)用于信號重構(gòu)的運行時間（CPU時間）對比，表中的數(shù)據(jù)是算法1與l1-ls，IST，GPSR_Basic，GPSR_Basic，MGPSR，運行 10次得到的平均值.數(shù)據(jù)顯示，算法 1與 l1-ls，IST，GPSR_Basic，GPSR_BB，MGPSR相比較，CPU時間最少，最快僅為IST的1/5.

圖1展示了原始信號，通過GPSR算法重構(gòu)的信號和算法1重構(gòu)的信號，其中的GPSR算法指的是單調(diào)的GPSR_BB算法.從圖形上可以看到：算法1得到的均方誤差（MSE） 較小，重構(gòu)的信號與原始信號接近，與單調(diào)的GPSR_BB算法重構(gòu)的信號相比，非零元素較少，信號的稀疏性較好，所以相比較GPSR_BB重構(gòu)的信號，算法1更有效.

圖1 從上到下依次為：原始信號，GPSR重構(gòu)的信號，算法1重構(gòu)的信號Fig.1 From top to bottom：original signal，reconstruction via the minimization of（2.1）obtained by GPSR，reconstruction via the minimization of（2.1） obtained by algorithm 1

圖2展示的是 GPSR_Basic，GPSR_BB，MGPSR，和算法1得到目標(biāo)函數(shù)值隨著迭代次數(shù)和CPU時間的變化.從圖形上看，目標(biāo)函數(shù)值隨著迭代次數(shù)和CPU時間的變化，算法1的迭代次數(shù)和CPU時間都相對較少.

圖3顯示了MSE隨著迭代次數(shù)和CPU時間的變化，從圖形上看，MSE隨著算法1的迭代次數(shù)和CPU時間的變化較快，算法1迭代次數(shù)和CPU時間都相對較少.

2.2 圖像復(fù)原試驗

圖像復(fù)原試圖利用退化現(xiàn)象的某種先驗知識復(fù)原被退化的圖像.復(fù)原技術(shù)是面向退化模型，并且采取相反的過程進行處理，以便恢復(fù)出原圖像.

下面將用算法1測試一些典型的測試圖像，這些圖像選自文獻[4]和文獻[6]，并與[6]中的算法結(jié)果進行了比較.在數(shù)值試驗中，通過離散小波變換得到稀疏圖像x，其中縮放因子和平移參數(shù)都選擇2j（j＞0的整數(shù)）的倍數(shù)，n=256，k=190，矩陣A∈Rk×n服從標(biāo)準(zhǔn)的高斯分布，且每行進行標(biāo)準(zhǔn)正交化，b=Ax，參數(shù)τ的選取為

圖2 目標(biāo)函數(shù)值隨迭代次數(shù)和CPU的變化情況Fig.2 The objective function plotted against iteration number and CPU time

圖3 MSE隨迭代次數(shù)和CPU的變化情況Fig.3 The MSE plotted against iteration number and CPU time

τ=0.02‖ATb‖∞

為了和MGPSR方法進行對比，終止條件與其終止條件保持一致，即

下面測試MGPSR和算法1這2種算法，測試圖像是256×256的圖像.

表2是MGPSR和算法1這2種算法圖像處理后的MSE，峰值信噪比（PSNR），CPU時間的對比.從 MSE，PSNR可以看出，用算法1復(fù)原的圖像與原始圖像最接近；從CPU時間可以看出，算法1速度要快，CPU時間最快大約僅為MGPSR的1/2.綜上所述，算法1明顯比MGPSR算法有效.

圖4對3個測試圖像進行了圖像處理，將算法1與MGPSR算法在圖像復(fù)原的質(zhì)量上進行了比較，可以看出，算法1復(fù)原的圖像的視覺效果明顯優(yōu)于MGPSR算法復(fù)原圖像的視覺效果.

表2 MGPSR和算法1兩種算法的比較Tab.2 The comparison of MGPSR and Algorithm1

3 結(jié)語

本文提出了一種新的應(yīng)用于稀疏信號重構(gòu)的加速動量梯度投影法，即把負梯度方向與動量項的凸組合作為搜索方向，步長采取LSD的步長選取規(guī)則.通過幾種算法數(shù)值試驗的比較，該方法在稀疏信號重構(gòu)的問題上明顯比其他幾種算法用時少，并且重構(gòu)的信號稀疏性較好；在圖像復(fù)原的問題上，復(fù)原的圖像明顯比其他算法復(fù)原的質(zhì)量好.本文的方法是有效的.但其收斂性和收斂速率還有待研究.

圖4 從左到右依次為：原始圖像，小波變換后的圖像，MGPSR復(fù)原的圖像，算法1復(fù)原的圖像Fig.4 From left to right：the original image，the image after wavelet transform，the image restored by MGPSR，the image restored by algorithm 1

[1]Figueiredo M A T，Nowak R D，Wright S J.Gradient projection for sparse reconstruction：application to compressed sensing and other inverse problems[C]//IEEE Journal of Selected Topics in Signal Processing.2007：586-597.

[2]Kim S J，Koh K，Lustig M，et al.An interior-point method for large-scale l1-regularized least squares[J].IEEE Journal of Selected Topics in Signal Processing，2007，1（4）：606-617.

[3]劉紫娟，李慧云，劉新為.外推系數(shù)帶參數(shù)的加速鄰近梯度算法[J].數(shù)值計算與計算機應(yīng)用，2016，37（3）：211-222.

[4]Zhou R，Niu L，Qi Z.Smoothing trust region for digital image restoration[C]//IEEE/WIC/ACM International Conference on Web Intelligence and Intelligent Agent Technology（WI-IAT）.IEEE，2015，3：40-43.

[5]Beck A，Teboulle M.A fast iterative shrinkage-thre holding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences，2009，2（1）：183-202.

[6]Ma G，Hu Y，Gao H.An accelerated momentum based gradient projection method for image deblurring[C]//IEEE International Conference on Signal Processing，Communications and Computing，2015.

[7]Figueiredo M A T，Nowak R D.A bound optimization approach to wavelet-based image deconvolution[C]//IEEE International Conference on Image Processing.2005：782-785.

[8]Amir B，Marc T.Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J].IEEE Transactions on Image Processing，2009，18（11）：2419-2434.

[9]Cao D D，Guo P.Blind image restoration based on wavelet analysis[C]//IEEE International Conference on Machine Learning and Cybernetics，2005，8：4977-4982.

[10]Rumelhart D E，Hinton G E，Williams R J.Learning internal representations by error propagation[R].California Univ San Diego La Jolla InstFor Cognitive Science，1985.

[11]畢亞倩，劉新為.求解界約束優(yōu)化的一種新的非單調(diào)譜投影梯度法[J].計算數(shù)學(xué)，2013，35（4）：419-430.

[12]Birgin E G，Martínez J M，Raydan M.Nonmonotone spectral projected gradient methods on convex sets[J].SIAM Journal on Optimization，2000，10（4）：1196-1211.

[13]Bertsekas D P.Projected Newton methods for optimization problems with simple constraints[J].SIAM Journal on control and Optimization，1982，20（2）：221-246.

[14]Facchinei F，Lucidi S，Palagi L.A truncated Newton algorithm for large scale box constrained optimization[J].SIAM Journal on Optimization，2002，12（4）：1100-1125.

[15]Doel K V D，Ascher U.The chaotic nature of faster gradient descent methods[J].Journal of Scientific Computing，2012，51（3）：1-22.

[責(zé)任編輯 楊 屹]

A new momentum gradient projection method for sparse reconstruction

WU Gaoyu1，LI Huiyun2
（1.School of Science,Hebei University of Technology,Tianjin 300401,China;2.School of Control Science and Engineering, Hebei University of Technology,Tianjin 300130,China）

A new momentum gradient projection method for sparse reconstruction is proposed by using the convex combination of the negative gradient direction and the momentum term as the search direction,and the same step length selection rule as the lagged steepest decent method.Different from the momentum based gradient projection method with fixed learning rate and momentum parameters,the proposed method employs dynamic selection of momentum parameters and step length,which accelerates its convergence.Experiment results demonstrate that the proposed method outperforms its competitors both in time efficiency and in the quality of signal reconstruction.

momentum gradient projection method;the lagged steepest descent method;l1regularization;signal reconstruction;image restoration

O224

A

1007-2373（2017）01-0059-06

10.14081/j.cnki.hgdxb.2017.01.010

2016-10-20

河北省自然科學(xué)基金（A2015202365）

武高玉（1990-），男，碩士研究生，shxwgy@sina.com.

：李慧云（1978-），女，講師，博士研究生，lhyhn@163.com.