馬 聰,劉 哲,甄小仙
西北工業(yè)大學(xué) 理學(xué)院,西安710129
信號(hào)采樣是聯(lián)系模擬信源和數(shù)字信息的橋梁。人們對(duì)信息的巨量需求造成了信號(hào)采樣、傳輸和存儲(chǔ)的巨大壓力。如何緩解這種壓力又能有效提取承載在信號(hào)中的有用信息是信號(hào)與信息處理中急需解決的問題之一。而近年來出現(xiàn)的壓縮感知理論(Compressed Sensing,CS)[1-2]提供了解決辦法,它突破了傳統(tǒng)的奈奎斯特采樣定律的限制,是一種可以同時(shí)完成采樣和壓縮的新型信號(hào)處理理論,其核心思想是:在某個(gè)變換域是稀疏的或者可壓縮的信號(hào),可以用一個(gè)與變換基不相關(guān)的觀測矩陣將變換所得高維信號(hào)投影到一個(gè)低維空間上,然后通過求解一個(gè)優(yōu)化問題從這些少量的投影中以高概率重構(gòu)出來。
在CS 理論中,信號(hào)的重構(gòu)基于范數(shù)優(yōu)化理論展開,數(shù)學(xué)意義上來講就是求解l0范數(shù)優(yōu)化問題,但該問題屬于非凸優(yōu)化問題[2],是NP 難的。文獻(xiàn)[3]證明了可以用l1范數(shù)代替l0范數(shù)求解得到同等的解,但該方法還需要較多的觀測樣本。文獻(xiàn)[4]構(gòu)造了lp(0 <p <1)范數(shù)實(shí)現(xiàn)算法,實(shí)驗(yàn)性地證明了基于lp(0 <p <1)范數(shù)的優(yōu)化算法在重構(gòu)信號(hào)效果及算法可靠性方面優(yōu)于l0范數(shù)優(yōu)化和l1范數(shù)優(yōu)化。文獻(xiàn)[5]從幾何學(xué)角度對(duì)不同范數(shù)的物理模型作出了解釋,并證明了lp范數(shù)在信號(hào)重構(gòu)方面的優(yōu)越性。
本文以非線性最優(yōu)化為方向研究lp范數(shù)優(yōu)化,提出了結(jié)合罰函數(shù)法和序列二次規(guī)劃方法(SQP)的lp范數(shù)優(yōu)化算法,實(shí)現(xiàn)了lp范數(shù)最優(yōu)化意義下CS 理論信號(hào)重構(gòu)。
設(shè)一個(gè)N 維有限長實(shí)信號(hào)x ∈RN,Φ 是一個(gè)M×N(M ?N)維的觀測矩陣,觀測值y=Φx 是一個(gè)M 維向量。在CS 理論框架下,原信號(hào)x 可以利用觀測向量y 通過求解如下問題得到重構(gòu):
其中‖·‖0表示向量中非零元素的個(gè)數(shù),也稱為l0范數(shù),但是該問題計(jì)算數(shù)值極不穩(wěn)定,很難利用現(xiàn)有的優(yōu)化方法有效求解。但是在滿足變換矩陣和觀測矩陣不相關(guān)的前提下,求解一個(gè)更加簡單的l1范數(shù)優(yōu)化問題代替l0范數(shù)優(yōu)化可以得到同等的解[3],即
這樣的改變將一個(gè)難以解決的非凸優(yōu)化問題轉(zhuǎn)化成了一個(gè)凸優(yōu)化問題,能夠方便地簡化為線性規(guī)劃進(jìn)行求解。但l1范數(shù)優(yōu)化還存在計(jì)算復(fù)雜度較高、需要觀測樣本較多和數(shù)據(jù)間有較大冗余性等缺陷。為了進(jìn)一步改進(jìn)重構(gòu)算法性能,直接最小化原信號(hào)的非凸lp(0 <p <1)范數(shù)[4],即
能夠很好地從數(shù)學(xué)原理上逼近信號(hào)重構(gòu)的原問題,不僅可以在一定程度上減少數(shù)據(jù)之間的冗余,而且大大減少準(zhǔn)確重構(gòu)原信號(hào)所需要的觀測數(shù)量。雖然lp范數(shù)優(yōu)化的研究成果很少,但是在信號(hào)重構(gòu)方面已經(jīng)表現(xiàn)出了巨大的優(yōu)勢,有效提高了基于CS 理論的信號(hào)重構(gòu)效果。
目前,最典型的lp(0 <p <1)范數(shù)實(shí)現(xiàn)算法是迭代加權(quán)最小二乘(IRLS)法[6],文獻(xiàn)[7]利用該算法從少量線性觀測中有效地重構(gòu)原矩陣信號(hào),并證明了任意由觀測向量唯一確定的稀疏信號(hào)均可以通過lp范數(shù)優(yōu)化方法精確重構(gòu)。在IRLS 算法基礎(chǔ)上,Rick Chartrand[8]提出了更為高效的參數(shù)規(guī)則化的IRLS 算法:
采用古典歐拉-拉格朗日乘子法求解式(4)得到其解:
Qn是一個(gè)對(duì)角矩陣,其元素為1/ωi。在迭代過程中,通過對(duì)參數(shù)ε 規(guī)則化處理,提高了對(duì)稀疏信號(hào)的重構(gòu)能力。
罰函數(shù)法[9]是一種求解非線性規(guī)劃問題的數(shù)值解法,它是將有約束的最優(yōu)化問題通過罰因子的選擇變?yōu)橐幌盗星罅P函數(shù)的無約束極小值問題。
以最小化lp范數(shù)為目標(biāo),文獻(xiàn)[10]從罰函數(shù)優(yōu)化角度提出了一種lp范數(shù)迭代重構(gòu)算法。令目標(biāo)函數(shù)為:f(x)=,約束條件轉(zhuǎn)化為:Φx-y=0,罰函數(shù)定義為:
通過增大罰因子,迫使迭代值進(jìn)入可行域并逼近最優(yōu)值。上式最優(yōu)解的迭代計(jì)算式為:
其中,mk為罰因子,Π(xk)=diag(|xk(i)p-2|) 。通過對(duì)目標(biāo)函數(shù)進(jìn)行非約束轉(zhuǎn)換和梯度修正,實(shí)現(xiàn)了對(duì)稀疏信號(hào)的精確重構(gòu)。但是隨著罰因子的增大,會(huì)導(dǎo)致約束條件變得很壞,從式(6)也可以看出,當(dāng)mk→∞時(shí),為了求出最優(yōu)解,必然要求Φx-y →0,這就給后期的計(jì)算帶來巨大的壓力,使得數(shù)值性能不好。因此在重構(gòu)稀疏度較高的信號(hào)時(shí),該方法不夠穩(wěn)定且計(jì)算復(fù)雜度較高。
序列二次規(guī)劃方法(Sequential Quadratic Programming,SQP)[11-12]是求解無約束最優(yōu)化問題的牛頓法和擬牛頓法對(duì)約束最優(yōu)化問題的推廣,其基本思想是每一步迭代通過求解一個(gè)二次規(guī)劃子問題來確定一個(gè)下降方向,以減少價(jià)值函數(shù)來取得步長,重復(fù)步驟直到求得原問題的解。
以非線性最優(yōu)化為方向?qū)崿F(xiàn)lp范數(shù)的優(yōu)化,采用SQP方法來求解lp范數(shù):
由問題式(9)取極小值的KT 條件[13]得到:
其中?f(x)=|p|diag(|x|p-2)x ,A(x)是約束函數(shù)(Φx-y)2的雅可比矩陣,采用牛頓法求解式(9),對(duì)于給定的點(diǎn)zk=(xk,μk),牛頓法迭代格式為:
這里pk=(dk,vk)滿足線性方程組:
其中τ >0 且充分小,式(12)等價(jià)為:
上式可以方便地利用二次規(guī)劃問題求解。
SQP 方法因具有很好的收斂性被廣泛應(yīng)用于求解約束最優(yōu)化問題,定理1[11]則從理論上闡明了該方法的收斂性。
定理1 假設(shè)目標(biāo)函數(shù)f(x) 與約束函數(shù)hj(x) 二階可微,問題式(8)的最優(yōu)解x*存在,雅可比矩陣A(x*)滿秩,存在μ*使得式(9)成立,如果初始點(diǎn)x0充分接近x*,則每次迭代二次規(guī)劃子問題的最優(yōu)解存在,迭代序列{xk}收斂于x*,且收斂速度是二次的。
罰函數(shù)法簡單直觀,隨著參數(shù)mk的增大,可以迫使初始點(diǎn)不斷向最優(yōu)解靠近,但是當(dāng)參數(shù)mk→∞時(shí)會(huì)使得迭代計(jì)算后期的數(shù)值性能不好。
考慮到以上兩種方法的特點(diǎn),提出結(jié)合罰函數(shù)與序列二次規(guī)劃的lp范數(shù)優(yōu)化算法(簡記為P-SQP 算法),首先利用罰函數(shù)法對(duì)初始點(diǎn)進(jìn)行預(yù)估計(jì),在這個(gè)過程中,不需要過多增大罰因子,僅是利用罰函數(shù)法做一個(gè)初步的計(jì)算,迫使初始點(diǎn)進(jìn)入可行域,然后采用SQP 方法求解計(jì)算。因?yàn)樯弦徊接?jì)算已經(jīng)保證了矩陣的正定性,根據(jù)定理1 利用SQP 方法二次收斂的優(yōu)點(diǎn)進(jìn)行后續(xù)計(jì)算可以快速求得原問題的最優(yōu)解。因此,P-SQP 算法結(jié)合了罰函數(shù)與SQP 兩種方法的優(yōu)點(diǎn)并克服了各自的不足,可以有效實(shí)現(xiàn)稀疏信號(hào)的重構(gòu)。
P-SQP 具體算法步驟如下:
步驟1選取初值x0∈Rn,μ0∈Rm,ρ,γ ∈(0,1),0 <ε ?1,誤差e,令k=0。
步驟2將x0代入式(7)求解,當(dāng)?shù)祒滿足‖Φx-y‖2<e 時(shí)停止計(jì)算,令x0=x。
步驟3計(jì)算‖ ?L(xk,μk) ‖,若‖ ?L(xk,μk) ‖≤ε ,停止計(jì)算,否則轉(zhuǎn)步驟4。
步驟4求解式(14),得dk和,并置vk=-μk-
步驟5令mk是使下式成立的最小非負(fù)整數(shù)m:
令αk=ρmk。
步驟6令xk+1=xk+αkdk,μk+1=μk+αkvk,k=k+1 返回步驟3。
下面通過實(shí)驗(yàn)將P-SQP 算法與IRLS 算法[8]和文獻(xiàn)[10]算法進(jìn)行比較分析,證明P-SQP 算法在一定程度上提升了對(duì)稀疏信號(hào)的重構(gòu)能力。
第一個(gè)實(shí)驗(yàn)是研究初值選取對(duì)P-SQP 算法收斂性的影響。實(shí)驗(yàn)信號(hào)x 選取長度為256 的稀疏信號(hào),觀測矩陣Φ 取100×256 的獨(dú)立同分布高斯矩陣,信號(hào)x 的稀疏度K依次為K=30,32,…,70,對(duì)每一個(gè)K 值,隨機(jī)生成100 個(gè)高斯隨機(jī)稀疏信號(hào)并采用P-SQP 算法做信號(hào)重構(gòu)實(shí)驗(yàn)(以下實(shí)驗(yàn)中,當(dāng)最優(yōu)解滿足‖-x‖2<10-3時(shí)視為精確重構(gòu))。算法第一步是用罰函數(shù)法對(duì)初始值進(jìn)行預(yù)估計(jì),使得迭代初值x0滿足‖ Φx0-y‖2<e,誤差e 的大小反映了x0和原信號(hào)x 的接近程度。圖1 給出了對(duì)于3 個(gè)不同大小的e 值,P-SQP 算法得出的精確重構(gòu)概率曲線。
圖1 不同e 值對(duì)應(yīng)的重構(gòu)概率圖
可以看出,e 值越小,即x0和x 越接近,P-SQP算法對(duì)信號(hào)的精確重構(gòu)概率越高,同時(shí)也證明了本文提出的P-SQP算法是有效可行的。
第二個(gè)實(shí)驗(yàn)是比較IRLS 算法,文獻(xiàn)[10]算法,P-SQP 算法對(duì)稀疏信號(hào)的精確重構(gòu)效率及平均計(jì)算時(shí)間,實(shí)驗(yàn)信號(hào)及觀測矩陣的選取同實(shí)驗(yàn)1,圖2 和圖3 分別給出了隨著信號(hào)稀疏度的增加,三種算法對(duì)信號(hào)x 的重構(gòu)概率曲線和平均計(jì)算時(shí)間曲線。
圖2 三種算法精確重構(gòu)信號(hào)概率圖
圖3 三種算法重構(gòu)信號(hào)平均計(jì)算時(shí)間曲線圖
圖4 P-SQP 算法重構(gòu)圖像
表1 三種算法重構(gòu)圖像結(jié)果
從圖2 和圖3 可以看出,對(duì)于相同的稀疏信號(hào),P-SQP算法的重構(gòu)概率是最高的;同時(shí),在信號(hào)稀疏度較低情況下,P-SQP 算法的收斂速度最快,計(jì)算時(shí)間是最少的。證明了P-SQP 算法充分結(jié)合了罰函數(shù)與SQP 兩種方法的優(yōu)點(diǎn),能夠有效地重構(gòu)稀疏信號(hào)。
第三個(gè)實(shí)驗(yàn)是研究P-SQP算法對(duì)稀疏圖像的重構(gòu)效果,首先對(duì)圖像施用水平梯度稀疏化算子得到稀疏圖像x1,采用P-SQP 算法對(duì)其進(jìn)行采樣重構(gòu)得到稀疏重構(gòu)圖像x2,然后利用水平梯度逆算子得到重構(gòu)圖像x。
選取大小為256×256 的Phantom 圖像和350×350 的MRI-Liver 圖像作為實(shí)驗(yàn)圖像,其水平梯度圖像的列最大稀疏度分別為K1=14,K2=148。采用P-SQP 方法對(duì)圖像進(jìn)行重構(gòu),實(shí)驗(yàn)結(jié)果如圖4 所示。
對(duì)于Phantom 圖像,當(dāng)采樣率分別為14.1%、15.6%時(shí),重構(gòu)圖像與原圖像的信噪比分別為66.341 0、72.019 6,均方誤差分別為0.015 2、0.004 1;而當(dāng)采樣率提高到17.2%時(shí),重構(gòu)圖像與原圖像的信噪比為88.764 5,均方誤差為8.642 4×10-5,達(dá)到了近似精確重構(gòu)。對(duì)稀疏度較大的MRILiver 圖像進(jìn)行重構(gòu),隨著采樣率的提高也得到了同樣的實(shí)驗(yàn)結(jié)果。證明了本文算法可以有效重構(gòu)稀疏圖像。
為比較IRLS 算法,文獻(xiàn)[10]算法,P-SQP 算法對(duì)稀疏圖像的重構(gòu)效果,分別取不同的觀測數(shù)進(jìn)行采樣,表1 給出了三種算法在不同采樣率下重構(gòu)圖像的客觀結(jié)果。
從表1 可以看出,當(dāng)達(dá)到一定的采樣率時(shí),對(duì)于不同的稀疏圖像,三種算法均可以實(shí)現(xiàn)精確重構(gòu)。在相同采樣率下,P-SQP 算法得到的重構(gòu)圖像與原圖像的信噪比更高,均方誤差也更小。并且隨著采樣率提高,可以實(shí)現(xiàn)對(duì)圖像的精確重構(gòu)時(shí),P-SQP 算法的計(jì)算時(shí)間是最少的。進(jìn)一步驗(yàn)證了本文算法的可行性與有效性。
lp范數(shù)問題是信號(hào)優(yōu)化重構(gòu)的一個(gè)新的方向,對(duì)它的研究將推動(dòng)壓縮感知理論在壓縮方面的應(yīng)用,但是lp范數(shù)屬于非凸函數(shù)優(yōu)化問題,目前的研究成果不是很多,很多數(shù)學(xué)問題還需要進(jìn)一步解決。本文沿著非線性最優(yōu)化方向,采用典型的約束優(yōu)化方法來求解lp(0 <p <1) 范數(shù)優(yōu)化問題。提出了結(jié)合罰函數(shù)和序列二次規(guī)劃方法的P-SQP 優(yōu)化算法,并進(jìn)行了實(shí)驗(yàn)分析,證明了該算法是一種有效的lp范數(shù)求解算法。如何選取初值或?qū)ζ漕A(yù)估計(jì)使其進(jìn)入可行域,使得該算法更具有應(yīng)用性還需進(jìn)一步研究。
[1] Donoho D.Compressed sensing[J].IEEE Trans on Information Theory,2006,52(4):1289-1306.
[2] Candès E.Compressive sampling[C]//International Congress of Mathematicians,Madrid,Spain,2006:1433-1452.
[3] Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Review,2001,43(1):129-159.[4] Chartrand R.Exact reconstruction of sparse signals via nonconvex minimization[J].IEEE Signal Process Lett,2007,14:707-710.
[5] Baraniuk R.A lecture on compressive sensing[J].IEEE Signal Processing Magazine,2007,24(4):118-121.
[6] Rao B D,Kreutz-Delgado K.An affine scaling methodology for best basis selection[J].IEEE Trans on Signal Process,1999,47:187-200.
[7] Fornasier M,Rauhut H,Ward R.Low-rank matrix recovery via iteratively reweighted least squares minimization[J].Communications on Pure and Applied Mathematics,2010,1:1-38.
[8] Chartrand R,Yin W.Iteratively reweighted algorithms for compressive sensing[C]//Proceedings of the 33rd International Conference on Acoustics,Speech,and Signal Processing(ICASSP),2008.
[9] 傅英定,成孝予,唐應(yīng)輝.最優(yōu)化理論與方法[M].北京:國防工業(yè)出版社,2008.
[10] 劉哲,楊揚(yáng).一種新的基于壓縮感知理論的稀疏信號(hào)重構(gòu)算法[J].光電子激光,2011,22(2).
[11] 孫文瑜,徐成賢,朱德通.最優(yōu)化方法[M].北京:高等教育出版社,2010.
[12] Gill P E,Wong E.Sequential quadratic programming methods,Technical Report NA-10-03[R].UCSD Department of Mathematics,2010.
[13] Jeyakumar V,Srisatkunarajah S.Geometric conditions for Kuhn-Tucker sufficiency of global optimality in mathematical programming[J].European Journal of Operational Research,2009,194:363-367.