壓縮感知中一種改進(jìn)內(nèi)點(diǎn)算法的研究

董 騰，楊 帆，潘國峰

(1.中國人民解放軍93704部隊(duì)，北京 101100；2.河北工業(yè)大學(xué)信息工程學(xué)院，天津 300401)

?

壓縮感知中一種改進(jìn)內(nèi)點(diǎn)算法的研究

董 騰1，2，楊 帆2，潘國峰2

(1.中國人民解放軍93704部隊(duì)，北京 101100；2.河北工業(yè)大學(xué)信息工程學(xué)院，天津 300401)

提出一種基于內(nèi)點(diǎn)法的改進(jìn)重構(gòu)算法，嘗試用專門的內(nèi)點(diǎn)算法解決稀疏重構(gòu)問題。首先，在內(nèi)點(diǎn)法基礎(chǔ)上引入預(yù)處理算子重新設(shè)計(jì)來避免牛頓方程系統(tǒng)的構(gòu)造，使矩陣擁有良好的可調(diào)性；其次，利用稀疏矩陣的矩陣特性簡(jiǎn)化矩陣矢量增量。仿真實(shí)驗(yàn)結(jié)果表明改進(jìn)的內(nèi)點(diǎn)算法對(duì)實(shí)際問題的處理是有效且優(yōu)于其他算法的。

壓縮感知；內(nèi)點(diǎn)法；稀疏矩陣；預(yù)處理；共軛梯度

0 引言

壓縮感知[1-3](Compressed Sensing，CS)理論是21世紀(jì)初提出的一種全新的信息獲取理論。該理論的提出能夠有效應(yīng)對(duì)目前大信息量數(shù)據(jù)傳輸和存儲(chǔ)成本極高的問題，因此理論一經(jīng)提出便引起國內(nèi)外相關(guān)學(xué)者的廣泛關(guān)注，基于壓縮感知的數(shù)據(jù)重構(gòu)理論、算法不斷被提出。

基于壓縮感知理論的重構(gòu)算法包括凸優(yōu)化算法、方向?qū)?yōu)算法和貪婪算法等[4]。凸優(yōu)化類算法通過凸優(yōu)化問題的極小化算法以極大概率逼近目標(biāo)函數(shù)以此重構(gòu)出原始信號(hào)，此類算法所需觀測(cè)次數(shù)少，重構(gòu)精度高，但計(jì)算復(fù)雜度高。目前主要有內(nèi)點(diǎn)法、基追蹤(BP)法和迭代閾值法等。

本文利用內(nèi)點(diǎn)算法[5]原始對(duì)偶特性，引入預(yù)處理共軛梯度對(duì)內(nèi)點(diǎn)法加以改進(jìn)，降低了算法復(fù)雜度,提高了算法的效率，在處理稀疏重構(gòu)尤其是大規(guī)模重構(gòu)問題時(shí)達(dá)到較好的重構(gòu)效果。目前，較先進(jìn)的重構(gòu)算法有定點(diǎn)連續(xù)有效集法(FPC-AS)[6]和光譜投影梯度(SPGL1)算法[7]，PDCO算法[8]是另外一種基于內(nèi)點(diǎn)算法的改進(jìn)算法，在文中利用這三種算法作為本文算法的比較算子，來比較改進(jìn)算法的計(jì)算效率。

1 壓縮感知的基本問題

壓縮感知基本問題為求不完備線性方程的解(m

(1)

研究表明，在某些場(chǎng)合中x稀疏解的精確重構(gòu)問題可以極大概率通過求解下列基追蹤[9]問題得到解決：

(2)

在實(shí)際應(yīng)用中，式(1)右邊經(jīng)常會(huì)有噪聲干擾，因此表示為

(3)

式中e為誤差，e∈Rm。

(4)

式中τ為調(diào)節(jié)稀疏度和噪聲誤差上限的正標(biāo)量。

式(4)就是著名的基追蹤降噪法。實(shí)際問題總是具有較大的規(guī)模，現(xiàn)成的辦法使得如單純形法或內(nèi)點(diǎn)法無法實(shí)現(xiàn)。然而，在壓縮感知問題中出現(xiàn)的矩陣A展示出了良好的特點(diǎn)，在優(yōu)化算法中可以加以利用?；诖吮疚奶岢鲆环N專門的方法來解決這個(gè)問題。

定義新變量u和ν∈Rn，則有

|xi|=ui+vi,?i=1,2,…,n

(5)

ui=max(xi，0)，νi=max(-xi，0)。則1范數(shù)形式變?yōu)?/p>

(6)

u，ν≥0，ln∈Rn成為一個(gè)列向量。用上述線性方法，解決式(4)中的約束光滑重構(gòu)：

(7)

式中：z=[u;ν]∈R2n；FT=[A,-A]∈Rm×2n。

一旦變量u和ν最優(yōu)解被找到，那么x的初始值就可通過計(jì)算下式得到：

x=u-v

線性化的意義在于初始化式(4)的維數(shù)是雙數(shù)，有2n個(gè)非負(fù)的約束條件加入。通過加入約束條件z≥0,使得線性搜索的每一步沿著負(fù)梯度方向并且新的迭代被投影到可行域。

2 內(nèi)點(diǎn)算法的改進(jìn)研究

2.1 內(nèi)點(diǎn)法的原始對(duì)偶問題

非光滑基追蹤式(2)和基追蹤去噪式(4)最優(yōu)化問題可以重新分別表示為等效線性和二次凸規(guī)劃問題。這些通過對(duì)目標(biāo)函數(shù)非光滑1范數(shù)規(guī)范化而得到。

對(duì)于原始問題式(7),它的對(duì)偶:

(8)

z,s≥0

在每一步原始對(duì)偶內(nèi)點(diǎn)法[5]應(yīng)用原始對(duì)偶對(duì)式(7)和式(8),相應(yīng)牛頓方向(Δz，Δs)的計(jì)算通過下邊系統(tǒng)的線性方程來解決:

(9)

式中：S和Z是s和z在對(duì)角線上的對(duì)角陣；I2n表示2n階的單位陣。

(10)

式中:μ=z′ts/(2n)是內(nèi)點(diǎn)法的約束項(xiàng);σ為定心參數(shù)，0≤σ≤1。

在稀疏矩陣結(jié)構(gòu)中，式(9)中的對(duì)偶變量Δs被除去，得到：

(Θ-1+FFT)Δz=fz+Z-1fs

(11)

Δs=Z-1fs-Θ-1Δz

(12)

式中Θ=S-1Z∈R2n×2n。

簡(jiǎn)化的牛頓系統(tǒng)式(11)又稱增廣系統(tǒng)，可由只有約束矩陣F存在的預(yù)處理迭代法處理。

2.2 預(yù)處理共軛梯度法解決原對(duì)偶問題

式(11)有一個(gè)對(duì)稱正定矩陣，在稀疏矩陣中可用共軛梯度(CG)法[10]來處理?？墒钱?dāng)矩陣是病態(tài)或特征值不集中時(shí)，CG法的收斂性很慢。在本節(jié)，討論一種有效的等量對(duì)角矩陣調(diào)節(jié)器來解式(11)的問題。

式(11)里提出的預(yù)調(diào)節(jié)器是利用式(11)CS矩陣一般特性和Θ矩陣的近似最優(yōu)解的通用屬性的。在式(7)、式(8)中原始對(duì)偶對(duì)的符號(hào)表示，變量s∈R2n是一個(gè)非負(fù)性約束z≥0的拉格朗日乘數(shù)。因此，在最優(yōu)解處sjzj=0(j=1，2，…，2n)。內(nèi)點(diǎn)法通過微擾這個(gè)約束條件sjzj=μ使收斂向著最優(yōu)解，μ是內(nèi)點(diǎn)法的約束項(xiàng)，使μ逐漸減少擾動(dòng)直至0。在最優(yōu)值j∈{1，2，…，2n}處，被分解為2部分：

(13)

這決定了約束條件的約束能力。這種分割對(duì)對(duì)角縮放矩陣Θ=S-1Z產(chǎn)生極大的不良影響：當(dāng)μ趨向于0時(shí)，指數(shù)j∈Β，Θj趨向無窮大，j∈Ν，Θj趨向于0。

前面提到z=[u;v]，u,v分別是向量x的正負(fù)部分。對(duì)于稀疏信號(hào)，只有k(k?2n)的非零部分是最優(yōu)解。u中將提供一個(gè)正的非零解，v將提供一個(gè)負(fù)的非零解。在Β中最優(yōu)解是k。因此，內(nèi)點(diǎn)法的后期迭代為

(14)

回到式(12)提出的預(yù)處理問題，它的矩陣為

H=Θ-1+FFT

(15)

矩陣Θ近似最優(yōu)解由式(14)表示。矩陣Θ-1有許多項(xiàng)，在內(nèi)點(diǎn)法達(dá)到最優(yōu)解[11]之前只有很少部分是完備解。設(shè)一個(gè)數(shù)C?1將Θ-1分解成不同部分：

(16)

l只是Θ-1里少數(shù)項(xiàng)的數(shù)目，與最優(yōu)解中的稀疏的k不同。在l

(17)

預(yù)調(diào)節(jié)器是基于ATA的近似，近似由封閉的擴(kuò)展單位陣ρIn得到，ρ=m/n：

(18)

為了簡(jiǎn)化對(duì)這個(gè)預(yù)調(diào)節(jié)器的分析，首先考慮n×n的矩陣H和P而不是式(17)、式(18)中所說的2n×2n的形式。下邊的引理建立了預(yù)處理矩陣P-1H在不分塊部分的頻譜特性：

引理1定義矩陣H：

H=Θ-1+ATA

式中:Θ=diag(Θ1,Θ2,…,Θn);Θ為n×n對(duì)角線陣Θj>0;A為m×n矩陣,m≤n/2。

P=Θ-1+ρIn，ρ=m/n

在1處收斂，也就是說：

(19)

如果矩陣A有近似正交行，可得：

描述式(11)中預(yù)處理矩陣P-1H的頻譜特性。使式(17)、式(18)中的H，P為分塊矩陣，則預(yù)處理矩陣P-1H有：n的多重特征值為1；其余n的特征值由引理1中Θ=Θu+Θv來定義。得到了P-1H在1附近收斂的特征值。因此，最終目的是希望應(yīng)用于式(11)的共軛梯度法的迭代在使用式(18)的預(yù)處理算子P時(shí)可以盡量少。在式(17)中H的矩陣A特征值λ(H)和λ(P-1H)的聚集是一個(gè)有標(biāo)準(zhǔn)正交行的離散余弦變換(DCT)矩陣，AAT=I。設(shè)數(shù)據(jù)規(guī)模為m=210，n=212，稀疏度k=51。得到如圖1所示特征值收斂分布。在圖1(a)中，顯示特征值λ(H)的聚集性。每一條垂線顯示特征值λ(H)在稀疏矩陣內(nèi)點(diǎn)法共軛梯度調(diào)用中的分布變化趨勢(shì)?？梢钥吹近c(diǎn)的聚集性隨著稀疏矩陣內(nèi)點(diǎn)法的優(yōu)化而變壞，而預(yù)處理矩陣P-1H特征值卻正相反，如圖1(b)。尤其是隨著稀疏矩陣內(nèi)點(diǎn)法的推進(jìn)特征值λ(P-1H)開始在1附近收斂[12-13]。

圖1 矩陣H和P-1H特征值的收斂性

2.3 改進(jìn)的內(nèi)點(diǎn)算法

迭代法(預(yù)處理共軛梯度法)處理式(11)、式(12)，在計(jì)算復(fù)合方向每一項(xiàng)時(shí)是近似相同的。為了避免過度計(jì)算，在每次迭代之前校正方向，稍稍偏離預(yù)測(cè)器指向中央的方向而只在必要時(shí)按照校正器指向方向。如同原對(duì)偶內(nèi)點(diǎn)法中長(zhǎng)變量一樣，這保證了每一次迭代目標(biāo)函數(shù)迅速減少而算法保持指向中心的路徑盡量小。當(dāng)許多有誤差的預(yù)測(cè)方向被執(zhí)行了，那么原對(duì)偶迭代開始向著可行域的邊界靠近。這導(dǎo)致了連續(xù)迭代的步長(zhǎng)變小。當(dāng)這種情況發(fā)生時(shí)，一個(gè)強(qiáng)有力的向心校正器開始發(fā)揮作用，使下一步迭代向著中心附近保證下一步迭代取值足夠大。理論上，原對(duì)偶方向步長(zhǎng)值要遠(yuǎn)離0以使全局收斂得到保證。這使得原對(duì)偶內(nèi)點(diǎn)法只需少量的迭代就能夠快速收斂。

改進(jìn)內(nèi)點(diǎn)算法(計(jì)為C-IPM)流程如下：

對(duì)于k=1，2，…有從zk到zk+1和從sk到sk+1的迭代：

While式(7)and式(8)對(duì)偶間隙≥εdo

σ=σ2

elseσ=σ1

end if

(預(yù)測(cè)器)

(校正器)

σ=σ3

end if

end while。

3 仿真實(shí)驗(yàn)

將通過仿真實(shí)驗(yàn)分析本文所提算法的性能，驗(yàn)證改進(jìn)算法在計(jì)算時(shí)間、重構(gòu)效率上與現(xiàn)有先進(jìn)算法的優(yōu)勢(shì)。所有實(shí)驗(yàn)都在MATLAB中運(yùn)行調(diào)試。都使用MATLAB R2010b 64-bit在CORE i7 下的win7系統(tǒng)中運(yùn)行。

表1是一維稀疏信號(hào)的重構(gòu)試驗(yàn)結(jié)果對(duì)比，選擇(n/2)×n形式的矩陣，信號(hào)長(zhǎng)度n變化從500到30 000，稀疏度k=51，通過不同算法得到的運(yùn)行時(shí)間的對(duì)比圖如表1所示。將本文所提改進(jìn)算法記為C-IPM法。實(shí)驗(yàn)表明，算法的運(yùn)行時(shí)間在矩陣規(guī)模在2 500×5 000以下時(shí)，本文所提改進(jìn)算法與現(xiàn)有先進(jìn)算法和基于內(nèi)點(diǎn)算法的其他算法相比，所用時(shí)間優(yōu)勢(shì)相當(dāng)明顯，對(duì)比從500～30 000各實(shí)驗(yàn)采集點(diǎn)可發(fā)現(xiàn)矩陣規(guī)模在5 000以內(nèi)時(shí)，改進(jìn)算法的計(jì)算效率至少提高1倍以上，當(dāng)矩陣規(guī)模逐漸增大，這種優(yōu)勢(shì)變得越來越明顯，計(jì)算效率相比其他算法有了較大優(yōu)勢(shì)，甚至PDCO算法已經(jīng)在現(xiàn)有硬件條件下無法計(jì)算出結(jié)果，說明了本算法在處理大規(guī)模重構(gòu)問題時(shí)效果非常理想，對(duì)于目前數(shù)據(jù)信息采集和傳輸規(guī)模越來越大的實(shí)際困難，利用改進(jìn)內(nèi)點(diǎn)算法處理大規(guī)模、超大規(guī)模重構(gòu)問題有很大的現(xiàn)實(shí)意義。

表1 不同算法計(jì)算效率對(duì)比

圖2是不同噪聲干擾下一維稀疏信號(hào)的重構(gòu)實(shí)驗(yàn)，信號(hào)長(zhǎng)度n=500，稀疏度k=51，分別加入噪聲均方差為0、0.01、0.03、0.05的高斯白噪聲對(duì)信號(hào)進(jìn)行測(cè)量，相應(yīng)的重構(gòu)結(jié)果如圖2(a)一圖2(d)所示。由圖2實(shí)驗(yàn)結(jié)果可以看出，在噪聲均方差低于0.05的不同程度噪聲干擾下，改進(jìn)內(nèi)點(diǎn)算法可精確重構(gòu)出原始信號(hào)。

圖3比較了本文所提的改進(jìn)算法與原始數(shù)據(jù)相變特性，參數(shù)n固定為1 000。測(cè)量值m變量范圍100～900。對(duì)每9個(gè)測(cè)量值m最優(yōu)表示的稀疏度k變化從1至m，對(duì)于每個(gè)k，都進(jìn)行100次實(shí)驗(yàn)。由圖3可以看出改進(jìn)內(nèi)點(diǎn)法實(shí)驗(yàn)相變特性與理論上的平均相變特性曲線相重疊，說明本文所提改進(jìn)算法在重構(gòu)質(zhì)量上達(dá)到了完全重構(gòu)。

4 結(jié)束語

本文設(shè)計(jì)并提出了一種低復(fù)雜度的稀疏矩陣原對(duì)偶內(nèi)點(diǎn)法來解決壓縮感知領(lǐng)域提出的l1正規(guī)化問題。在原始對(duì)偶內(nèi)點(diǎn)法每次迭代中，迭代方向由共軛梯度法解線性系統(tǒng)得到。然而矩陣Θ-1+FFT作為算法收斂是一個(gè)病態(tài)的，因此，共軛梯度法計(jì)算緩慢。為了解決這一病態(tài)問題，本文為共軛梯度法設(shè)計(jì)出一個(gè)低成本的預(yù)調(diào)節(jié)器。數(shù)值實(shí)驗(yàn)證實(shí)了改進(jìn)算法在計(jì)算效率上有了明顯的提升，為解決大規(guī)模信息處理問題提供了一種新的嘗試。

(a)稀疏隨機(jī)陣(噪聲均方差：0)

(b)稀疏隨機(jī)陣(噪聲均方差：0.01)

(c)稀疏隨機(jī)陣(噪聲均方差：0.03)

(d)稀疏隨機(jī)陣(噪聲均方差：0.05)圖2 不同噪聲干擾下一維稀疏信號(hào)的重鉤

圖3 原始數(shù)據(jù)與改進(jìn)算法相變特性比較

[1] 張波，劉郁林，王開.稀疏隨機(jī)有限等距性質(zhì)分析.電子與信息學(xué)報(bào)，2014，36(1)：169-174.

[2] 寧方立，何碧靜，韋娟.基于lp范數(shù)的壓縮感知圖像重建算法研究. 物理學(xué)報(bào)，2013(17)：174-212.

[3] SUN H, YUAN H M, LU S M. A design of three-phase digital watt-hour power meter on SOPC platform.Proc. of International Conference on Information Technology and Computer Science：IEEE Press, 2009.

[4] 李珅，馬彩文，李艷，等.壓縮感知重構(gòu)算法綜述.紅外與激光工程，2013，42(S1)：225-232.

[5] 雍龍泉.線性規(guī)劃的原-對(duì)偶內(nèi)點(diǎn)算法數(shù)值實(shí)驗(yàn)初步.科學(xué)技術(shù)與工程，2007，7(18)：4576-4579.

[6] WEN Z, YIN W, GOLDFARB D, et al. A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization and continuation.Sci. Comput.32(4):1809-1831, 2010.

[7] BERG E V, FRIEDLANDER M P. Probing the Pareto frontier for basis pursuit solutions. SIAM Journal of Scientic Computing,2008,31(2):890-912.

[8] SAUNDERS M, KIM B. PDCO：Primal-dual interior method for convex objectives.Technical Report, Stanford University, 2002. http://www.stanford.edu/group/SOL/software/pdco.html.

[9] CHEN S, DONOHO D L, SAUNDERS M A. Atomic decomposition by basis pursuit.SIAM Journal on Scientic Computing, 20(1):33-61, 1998.

[10] 董曉亮，李郴良，唐清肝.一類Wolfe搜索下的共軛梯度法及其全局收斂性.廣西科學(xué)，2007，14(1)：44-46.

[11] 姜志俠，李軍，張珊.一個(gè)改進(jìn)的原始對(duì)偶內(nèi)點(diǎn)法.吉林大學(xué)學(xué)報(bào)，2009，47(4)：677-682.

[12] HAUPT J, BAJWA W, RAZ G, et al. Toepitz compressed sensing matrices with applications to sparse channel estimation.IEEE Transactions Information Theory, 2010, 56(11):5862-5875.

[13] LIU Y L, WANG K,HE J W. Signal recovery by compressed sensing in IR-UWB systems. Chinese Journal of Electronics, 2012, 21(2)：339-344.

Algorithm for Modified Interior Point Methods on Compressed Sensing

DONG Teng1,2,YANG Fan2,PAN Guo-feng2

(1.The 93704 Unit of the Chinese People’s Liberation Army,Beijing 101100,China;2.School of Information Engineering,Hebei University of Technology,Tianjin 300401,China)

To use a special interior point method solve the sparse reconstruction problem,a modified reconstructed method based on interior point method was proposed.First,to avoid construction of Newton equations system,a preconditioning operator was lead up which was based on the interior point method to make the matrix have good adjustability.Second,the characteristics of sparse matrix were used to simplify the matrix-vector incrementing.The experimental results show that the modified interior point method is effective and better than any other algorithms for practical problems.

compressed sensing;interior point method;sparse matrix;precondition;conjugate gradient

國家科技重大專項(xiàng)資助項(xiàng)目(2009ZX02308-004)

2014-07-30 收修改稿日期：2015-03-28

TP391

A

1002-1841(2015)06-0138-05

董騰(1987—)，碩士研究生，助理工程師，主要從事圖像處理、智能信息處理及信號(hào)分析等方面的研究。 E-mail:dongteng2006@sina.com 楊帆(1964—)，博士，教授，博士生導(dǎo)師，主要從事計(jì)算機(jī)視覺檢測(cè)技術(shù)、生物特征識(shí)別技術(shù)和圖像處理與模式識(shí)別方面的研究。