董 騰，楊 帆，潘國峰

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

?

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

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

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

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

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

0 引言

壓縮感知[1-3](Compressed Sensing，CS)理論是21世紀(jì)初提出的一種全新的信息獲取理論。該理論的提出能夠有效應(yīng)對目前大信息量數(shù)據(jù)傳輸和存儲成本極高的問題，因此理論一經(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)出原始信號，此類算法所需觀測次數(shù)少，重構(gòu)精度高，但計算復(fù)雜度高。目前主要有內(nèi)點法、基追蹤(BP)法和迭代閾值法等。

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

1 壓縮感知的基本問題

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

(1)

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

(2)

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

(3)

式中e為誤差，e∈Rm。

(4)

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

式(4)就是著名的基追蹤降噪法。實際問題總是具有較大的規(guī)模，現(xiàn)成的辦法使得如單純形法或內(nèi)點法無法實現(xiàn)。然而，在壓縮感知問題中出現(xiàn)的矩陣A展示出了良好的特點，在優(yōu)化算法中可以加以利用。基于此本文提出一種專門的方法來解決這個問題。

定義新變量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成為一個列向量。用上述線性方法，解決式(4)中的約束光滑重構(gòu)：

(7)

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

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

x=u-v

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

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

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

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

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

(8)

z,s≥0

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

(9)

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

(10)

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

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

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

(11)

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

(12)

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

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

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

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

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

(13)

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

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

(14)

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

H=Θ-1+FFT

(15)

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

(16)

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

(17)

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

(18)

為了簡化對這個預(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對角線陣Θ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時可以盡量少。在式(17)中H的矩陣A特征值λ(H)和λ(P-1H)的聚集是一個有標(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ào)用中的分布變化趨勢。可以看到點的聚集性隨著稀疏矩陣內(nèi)點法的優(yōu)化而變壞，而預(yù)處理矩陣P-1H特征值卻正相反，如圖1(b)。尤其是隨著稀疏矩陣內(nèi)點法的推進(jìn)特征值λ(P-1H)開始在1附近收斂[12-13]。

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

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

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

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

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

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

σ=σ2

elseσ=σ1

end if

(預(yù)測器)

(校正器)

σ=σ3

end if

end while。

3 仿真實驗

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

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

表1 不同算法計算效率對比

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

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

4 結(jié)束語

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

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

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

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

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

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

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

[2] 寧方立，何碧靜，韋娟.基于lp范數(shù)的壓縮感知圖像重建算法研究. 物理學(xué)報，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ī)劃的原-對偶內(nèi)點算法數(shù)值實驗初步.科學(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] 姜志俠，李軍，張珊.一個改進(jìn)的原始對偶內(nèi)點法.吉林大學(xué)學(xué)報，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

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

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

TP391

A

1002-1841(2015)06-0138-05

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