基于回溯的共軛梯度迭代硬閾值重構(gòu)算法

張雁峰，范西岸，尹志益，蔣鐵鋼

(1.廣東工業(yè)大學(xué) 機(jī)電工程學(xué)院，廣州 510006； 2.廣東科技學(xué)院 機(jī)電工程學(xué)院, 廣東 東莞 523083)(*通信作者電子郵箱1648421565@qq.com)

0 引言

在信號處理領(lǐng)域，傳統(tǒng)的奈奎斯特采樣方法的采樣頻率須高于實(shí)際信號最高頻率的2倍，如此高密度的采樣方式給數(shù)據(jù)量日益增加的信號處理過程增加了壓力。壓縮感知(Compressed Sensing, CS)，是由Donoho[1]提出來的一種亞采樣信息獲取理論。跟傳統(tǒng)的奈奎斯特采樣理論不同，CS采樣頻率可以遠(yuǎn)小于實(shí)際信號最高頻率的兩倍，即只需少量采樣便可精確地重構(gòu)原信號，從而大大降低了信號的存儲和運(yùn)輸?shù)拇鷥r(jià)。目前壓縮傳感理論已經(jīng)在信息論、醫(yī)療成像、模式識別、雷達(dá)探測、地質(zhì)勘探、圖像壓縮、圖像超分辨率重建等領(lǐng)域受到高度關(guān)注和應(yīng)用[2]。信號重構(gòu)算法關(guān)系著信號重構(gòu)速度和質(zhì)量，一直是CS的研究熱點(diǎn)。

對于長度為N的信號f，可以表示為一組正交基向量Ψ={ψi|i=1,2,…，N}的線性組合，即:

f=Ψx

(1)

其中:Ψ∈RN×N稱為稀疏基;x∈RN×1為f在Ψ變換域中的系數(shù)向量，如果x中只有s(s?N)個元素不為零(或遠(yuǎn)大于零，而其他元素接近于零)，則稱f是s-稀疏的。

用一個M×N(M

y=ΦΨx=Ax

(2)

其中A∈RM×N稱為傳感矩陣。文獻(xiàn)[3]證明了當(dāng)傳感矩陣A滿足s階有限等距性質(zhì)(Restricted Isometry Property, RIP)時(shí)，可以通過求解式(3)以很高的概率精確求得x：

(3)

s. t.y=Ax

求解式(3)的過程稱為CS信號重構(gòu)。貪婪算法作為一種CS重構(gòu)算法，由于具有操作簡單、計(jì)算復(fù)雜度低的特點(diǎn)而被廣泛應(yīng)用。正交匹配追蹤(Orthogonal Matching Pursuit, OMP)算法是一種典型的貪婪算法[4]，OMP算法每次迭代挑選一個原子，當(dāng)采樣數(shù)量較大或稀疏度較大時(shí)，算法復(fù)雜度急劇增加。文獻(xiàn)[5]使用共軛梯度法代替OMP算法中每次迭代所需的最小二乘運(yùn)算，在稀疏度較大的情況下優(yōu)勢明顯。文獻(xiàn)[6-8]在OMP算法的基礎(chǔ)上，每次迭代挑選多個原子，一定程度上減少了找到正確支撐集的時(shí)間，但是需要選擇合適的參數(shù)。文獻(xiàn)[9]提出的迭代硬閾值(Iterative Hard Thresholding, IHT)算法是一種基于l0范數(shù)的凸優(yōu)化算法，算法復(fù)雜度低，但是容易受傳感矩陣縮放影響導(dǎo)致不穩(wěn)定。文獻(xiàn)[10]使用利普希茨收斂條件將IHT算法中的固定步長轉(zhuǎn)化為一個下降序列，提出正則化 IHT(Normalized IHT, NIHT)算法，算法復(fù)雜度稍有增加，但是迭代次數(shù)大幅減少。IHT/NIHT的缺點(diǎn)是有可能導(dǎo)致某些支撐被重復(fù)選擇，從而導(dǎo)致其找到正確的支撐集需要較多的迭代次數(shù)[11]。文獻(xiàn)[12]提出共軛梯度 IHT(Conjugate Gradient IHT, CGIHT)算法, 該算法使用共軛梯度法逼近最優(yōu)解，減少了迭代次數(shù)，但是在實(shí)際過程中，由于誤差的累積和傳感矩陣的列原子不完全正交，導(dǎo)致前后迭代中尋優(yōu)方向不共軛，當(dāng)?shù)阶顑?yōu)解附近時(shí)還是可能出現(xiàn)類似梯度下降法的振蕩現(xiàn)象，導(dǎo)致算法收斂速度變慢。

文獻(xiàn)[13-15]中算法引入了回溯的思想，降低原子被重復(fù)選擇的概率，計(jì)算復(fù)雜度相對較低且重構(gòu)精度高。其中文獻(xiàn)[13]中的子空間追蹤(Subspace Pursuit, SP)算法和文獻(xiàn)[14]中的基于回溯的IHT算法(Backtracking-based IHT algorithm, BIHT)在每次迭代中取當(dāng)前迭代的s個支撐與前一次迭代的s個支撐組成大小為2s的候選集，并將觀測向量重新投影到候選集對應(yīng)的矩陣列所張成的空間，從投影值中選出s個支撐作為新的支撐集。文獻(xiàn)[15]算與SP算法/BIHT采用類似的回溯過程，不同的是其每次迭代從大小為3s的候選集中選出s個元素作為支撐集。SP算法和文獻(xiàn)[15]算法都是以O(shè)MP算法為基礎(chǔ)，且都采用回溯思想，兩種算法的重構(gòu)性能相當(dāng)。BIHT以IHT算法為基礎(chǔ)，通過定點(diǎn)迭代法求解凸優(yōu)化問題式(4)來求稀疏解，與SP算法和文獻(xiàn)[15]算法采用最小二乘求解相比，BIHT迭代次數(shù)更多。

s. t. ‖x‖0≤s

(4)

本文以BIHT為基礎(chǔ)，保留其每次迭代過程中的回溯過程，用梯度下降法和共軛梯度法結(jié)合的方式代替BIHT中單純使用梯度下降法作為尋優(yōu)方法，提出了一種基于回溯的共軛梯度迭代硬閾值算法(Backtracking-based Conjugate Gradient Iterative Hard Thresholding algorithm, BCGIHT)。通過一維信號和二維圖像的重構(gòu)仿真實(shí)驗(yàn)，與NIHT算法、BIHT、CGIHT算法、文獻(xiàn)[15]的算法作對比，驗(yàn)證了BCGIHT的性能。

1 基于回溯的共軛梯度迭代硬閾值算法

1.1 算法描述

BCGIHT主要策略是：一方面，結(jié)合BIHT的回溯思想，使得迭代過程中支撐集被快速準(zhǔn)確找到，確保算法重構(gòu)概率和重構(gòu)精度；另一方面，采用梯度下降法和共軛梯度法相結(jié)合的方式代替迭代硬閾值類算法所使用的單純的梯度下降法，進(jìn)一步提高算法的收斂速度。

首先在每次迭代中采用非線性算子Hs(α)獲得系數(shù)向量xn中絕對值最大的s個支撐，并取前一步迭代的支撐集與當(dāng)前s個支撐的并集作為候選集，將觀測向量y在候選集對應(yīng)的矩陣列所張成的空間重新做一次正交投影得到中間系數(shù)向量an，再篩選出an中絕對值最大的s個支撐作為新的支撐集。由于每次迭代中回溯了前一步迭代中xn-1的支撐，減少了某些支撐在第n-1次迭代中被設(shè)為0而在第n次迭代中重新被選擇的次數(shù)，使得支撐被反復(fù)選擇的概率大大降低，優(yōu)化了支撐集的選擇，從而可以快速找到正確支撐集。

優(yōu)化方法中梯度下降法具有線性收斂速度，共軛梯度法具有超線性收斂速度。由于誤差的累積及實(shí)際過程中傳感矩陣的列原子不完全正交，導(dǎo)致在迭代過程中，前后尋優(yōu)方向并不共軛，因此使用單純的共軛梯度法作為尋優(yōu)方法效果不好。BCGIHT在BIHT單純使用梯度下降法的基礎(chǔ)上引入共軛梯度法加速算法的收斂。以前后迭代所得支撐集是否相等作為準(zhǔn)則來確定尋優(yōu)方向：若前后支撐集相等，即Γn=Γn-1，說明由前后支撐集所對應(yīng)的觀測矩陣的列組成的矩陣相同，即AΓn=AΓn-1，因此AΓnTAΓn=AΓn-1TAΓn-1，這時(shí)可以繼續(xù)用共軛梯度方向進(jìn)行下一步尋優(yōu)；而當(dāng)前后迭代中的支撐集不相等時(shí)，上述兩個對稱正定矩陣AΓnTAΓn和AΓn-1TAΓn-1不再相等，前后尋優(yōu)方向dn和dn-1不再是共軛的，因此調(diào)整尋優(yōu)方向?yàn)樘荻认陆捣较?。相較于單純地使用梯度方向或者共軛方向的方法，尋優(yōu)方向更替策略能在每次迭代中自動調(diào)整尋優(yōu)方向，從而加速算法收斂。

1.2 算法流程

以下流程中:Hs(α)是一個非線性算子，表示取α中絕對值最大的s個元素;μ為迭代步長;Γn表示第n次迭代的索引集合;〈·,·〉表示求內(nèi)積;supp(·)表示向量非零元素索引;ε表示迭代停止誤差閾值。BCGIHT算法具體流程如下：

步驟1 初始化x0=0，殘差r0=y，初始梯度g0=ATy，初始尋優(yōu)方向d0=g0，初始支撐集Γ0=supp(Hs(g0))，μ0=〈g0,g0〉/〈Ag0,Ag0〉，x1=Hs(x0+μ0d0)，Γ1=supp (x1)。

步驟2n≥1時(shí)，更新梯度方向gn=AT(y-Axn)，判斷是否Γn=Γn-1，如果是，轉(zhuǎn)步驟4，否則，轉(zhuǎn)步驟3。

步驟3 梯度下降法。計(jì)算步長μn=〈gn Γ n,gn Γ n〉/〈AΓngnΓ n,AΓ ngnΓn〉，計(jì)算an=Hs(xn+μngn)。

步驟4 共軛梯度法。計(jì)算方向因子βn=〈gn Γ n,gnΓ n〉/〈g(n-1)Γ n-1,g(n-1)Γ n-1〉,更新dn=gn+βndn-1、μn=〈dnΓ n,gnΓ n〉/〈AΓ ndnΓ n,AΓndnΓn〉，計(jì)算an=Hs(xn+μndn)。

步驟6 更新殘差rn+1=y-Axn+1。

步驟7 判斷是否‖rn+1‖2<ε，若是，迭代停止；否則，n=n+1，返回步驟2繼續(xù)迭代。

2 算法性能分析

1)算法收斂性分析。BCGIHT采用梯度與共軛梯度交替的方式作為尋優(yōu)方向，由IHT算法和共軛梯度法的收斂性證明，可知在每次迭代中，殘差‖rn‖2<‖rn-1‖2，因此BCGIHT算法至少可以收斂到局部最優(yōu)解。

表1 不同重構(gòu)算法所需的迭代次數(shù)Tab. 1 Iterative number required for different algorithms

3 仿真結(jié)果與分析

為了驗(yàn)證BCGIHT的性能，本文使用一維信號和圖像進(jìn)行一系列仿真實(shí)驗(yàn)，并將該算法與具有代表性的NIHT算法、BIHT、CGIHT算法、文獻(xiàn)[15]算法作仿真對比。所有的實(shí)驗(yàn)都是在Intel i5-7500、8 GB DDR4內(nèi)存和Matlab 2014b上完成。

3.1 一維信號重構(gòu)

本次實(shí)驗(yàn)使用的測試信號長度為N=256，使用M×N的隨機(jī)高斯矩陣作為采樣矩陣，采樣數(shù)量M=128，每次實(shí)驗(yàn)重復(fù)做1 000次，以平均值作為實(shí)驗(yàn)結(jié)果。重構(gòu)成功率是指對同一實(shí)驗(yàn)進(jìn)行多次重復(fù)實(shí)驗(yàn)，其重構(gòu)均方根誤差在一定的誤差閾值λ內(nèi)的概率，本次實(shí)驗(yàn)將誤差閾值設(shè)為λ=E-12。在不同稀疏度下，NIHT、CGIHT、BIHT、BCGIHT和文獻(xiàn)[15]算法的重構(gòu)成功率對比結(jié)果如圖1所示。

從圖1可以看出,當(dāng)采樣次數(shù)固定，稀疏度越來越大時(shí)，所有算法的重構(gòu)成功率都會變小，這是由于稀疏度s變大，導(dǎo)致由s個大系數(shù)組成的支撐的不確定性增加導(dǎo)致的;NIHT算法在稀疏度比較小時(shí)，重構(gòu)成功率能保持較高水平，但是當(dāng)s>35時(shí)，重構(gòu)成功率急劇下降;BIHT和CGIHT相較于NIHT稍好，BCGIHT相對于BIHT和CGIHT有較大提高，在s=45時(shí)還能保持較高的重構(gòu)成功率;文獻(xiàn)[15]算法在s>40時(shí)重構(gòu)成功率急劇下降，這是由于稀疏度較大時(shí)，候選支撐集元素之間相互影響導(dǎo)致。

圖1 不同稀疏度下不同算法重構(gòu)成功率對比Fig. 1 Reconstruction success rate comparison of different algorithms with different sparsity

一維信號在不同稀疏度情況下不同算法重構(gòu)時(shí)間對比如圖2所示。由圖2可知，隨著稀疏度的增加，五種算法的重構(gòu)時(shí)間都趨于增加，且文獻(xiàn)[15]算法重構(gòu)時(shí)間增加得越來越明顯；BCGIHT重構(gòu)時(shí)間少于BIHT和文獻(xiàn)[15]算法25%以上,BCGIHT重構(gòu)時(shí)間介于BIHT和CGIHT算法之間;CGIHT算法所用時(shí)間最少，這是因?yàn)楸緦?shí)驗(yàn)設(shè)置的一維信號是嚴(yán)格稀疏信號，稀疏基為單位矩陣，不存在原子之間的相互干擾，AΓnTAΓn為嚴(yán)格對稱正定矩陣，因此CGIHT算法能以很少的迭代次數(shù)達(dá)到迭代停止條件。

圖2 不同稀疏度下不同算法重構(gòu)時(shí)間對比Fig. 2 Reconstruction time comparison of different algorithms with different sparsity

3.2 圖像重構(gòu)

本次實(shí)驗(yàn)采用常用的測試圖像Pepper(像素為512×512)驗(yàn)證BCGIHT重構(gòu)效果。首先，對原始圖像進(jìn)行小波變換，然后利用高斯隨機(jī)矩陣作為觀測矩陣進(jìn)行壓縮采樣，最后使用不同的算法重構(gòu)圖像，并比較分析重構(gòu)結(jié)果。

在采樣率為0.5(即M=0.5×512=256)、稀疏度為s=M/4=64時(shí),不同算法重構(gòu)出的Pepper圖像的直觀效果如圖3所示。其中圖3(b)～(f)的峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)分別為33.05 dB、36.34 dB、31.58 dB、36.15 dB、34.58 dB。

圖3 不同算法的重構(gòu)結(jié)果Fig. 3 Reconstruction results of different algorithms

為了比較在不同的采樣率下不同算法重構(gòu)質(zhì)量和重構(gòu)時(shí)間，實(shí)驗(yàn)設(shè)置了三組不同的采樣率，分別為0.3、0.5和0.7。用PSNR反映重構(gòu)質(zhì)量，其計(jì)算式為：

PSNR=10lg((N-1)2/MSE);

(5)

每種采樣率下不同的算法分別進(jìn)行100次蒙特卡洛實(shí)驗(yàn)，計(jì)算平均值作為結(jié)果。不同采樣率下各重構(gòu)算法的重構(gòu)PSNR和重構(gòu)時(shí)間對比如表2所示。從表2中可知，在采樣率為0.3和0.5時(shí)，BIHT和BCGIHT的PSNR略高于文獻(xiàn)[15]算法;在采樣率為0.7時(shí)，三種算法PSNR相當(dāng),相較而言，NIHT算法和CGIHT算法PSNR較小。在采樣率小于等于0.5時(shí)，BCGIHT重構(gòu)時(shí)間最短；CGIHT算法和BCGIHT的重構(gòu)時(shí)間均小于BIHT，這是因?yàn)锽IHT中有最小二乘計(jì)算，因此重構(gòu)時(shí)間相對較長，而BCGIHT雖然也需最小二乘計(jì)算，但是因?yàn)閷?yōu)方向改進(jìn)了，因此其迭代次數(shù)減少了，所以其重構(gòu)時(shí)間不僅沒有增加，反而減少了，且相較于BIHT，對于不同的采樣率，BCGIHT重構(gòu)時(shí)間減少50%以上；文獻(xiàn)[15]算法重構(gòu)時(shí)間介于BIHT和BCGIHT之間，符合理論分析。

在原圖上加上不同方差的高斯白噪聲時(shí)，不同算法重構(gòu)帶噪圖像的重構(gòu)精度對比如表3表示。對比表2和表3可知，相比于原圖重構(gòu)，CGIHT算法PSNR減少幅度較小，且噪聲方差越大時(shí)表現(xiàn)越明顯，這是因?yàn)镃GIHT算法沒有使用回溯方法挑選支撐集，不存在候選集原子之間相互干擾的情況;NIHT算法、BIHT、BCGIHT和文獻(xiàn)[15]算法抗噪性能相當(dāng);噪聲方差較大時(shí)，BCGIHT抗噪性能略好。

表2 不同采樣率下不同算法的重構(gòu)PSNR和重構(gòu)時(shí)間Tab. 2 Reconstruction PSNR and time for different algorithms at different sample rates

表3 采樣率為0.5時(shí)噪聲環(huán)境下不同算法重構(gòu)PSNRTab. 3 Reconstruction PSNR of different algorithms with Gaussian noise at the sample rate of 0.5

4 結(jié)語

本文在壓縮感知BIHT的研究基礎(chǔ)上，針對BIHT迭代次數(shù)多且重構(gòu)時(shí)間長的問題，提出BCGIHT并對其進(jìn)行了收斂性和算法復(fù)雜度分析:一方面，用回溯思想確保正確支撐集被快速找到;另一方面，采用梯度與共軛梯度結(jié)合使用的策略來逼近最優(yōu)解。一維信號和二維圖像的仿真實(shí)驗(yàn)結(jié)果表明BCGIHT在同精度條件下與BIHT及類似算法相比，迭代次數(shù)更少，重構(gòu)時(shí)間遠(yuǎn)低于BIHT，兼顧了算法重構(gòu)精度和重構(gòu)效率,且在不同強(qiáng)度噪聲環(huán)境下的重構(gòu)性能與BIHT及同類算法相當(dāng)，具有較好的抗噪性。