一類(lèi)信號(hào)恢復(fù)問(wèn)題的投影迭代硬閾值算法?

李慶龍， 曾雪迎

(中國(guó)海洋大學(xué)數(shù)學(xué)科學(xué)學(xué)院， 山東 青島 266100)

稀疏恢復(fù)在信號(hào)和圖像處理、機(jī)器學(xué)習(xí)等許多領(lǐng)域中都得到了廣泛的關(guān)注?；诓糠指道锶~系數(shù)恢復(fù)原始信號(hào)是稀疏恢復(fù)中的一類(lèi)重要問(wèn)題，在壓縮感知、磁共振成像等領(lǐng)域中具有重要的應(yīng)用價(jià)值[1]。該問(wèn)題需要尋找一個(gè)可以在傅里葉變換域內(nèi)擬合采樣不足的傅里葉系數(shù)的信號(hào)。稀疏恢復(fù)是典型的不適定反問(wèn)題。數(shù)據(jù)可稀疏逼近使得從部分傅里葉系數(shù)中恢復(fù)原始信號(hào)成為可能[2]。該問(wèn)題可以建模為L(zhǎng)0模最優(yōu)化問(wèn)題，但它是一個(gè)非凸不連續(xù)的NP難問(wèn)題，模型的理論分析和數(shù)值求解均非常困難。在不增加測(cè)量數(shù)據(jù)的情形下，已有文獻(xiàn)的解決方案通常是將L0模松弛為L(zhǎng)1模[2-4]或其它非凸但連續(xù)的稀疏性度量函數(shù)[5]。

雖然L1模對(duì)應(yīng)的凸優(yōu)化問(wèn)題具有快速算法[6-7],但其會(huì)在恢復(fù)信號(hào)中引入偏差，因此通常不能取得令人滿(mǎn)意的恢復(fù)效果。對(duì)于低維問(wèn)題，增加觀測(cè)傅里葉系數(shù)的數(shù)量可以帶來(lái)更準(zhǔn)確的解決方案，但它耗時(shí)且耗力，而且高維問(wèn)題將會(huì)挑戰(zhàn)設(shè)備硬件采樣和存儲(chǔ)的時(shí)間和物理成本，特別是在諸如磁共振成像等具體應(yīng)用背景中代價(jià)較大。將L0模松弛為非凸連續(xù)函數(shù)的方法雖然取得了很好的效果，但其適用范圍相對(duì)有限，目前主要用于求解稀疏信號(hào)的變換域系數(shù)，一般不能在信號(hào)域直接解決問(wèn)題。

已有文獻(xiàn)在構(gòu)造數(shù)據(jù)稀疏變換方面做了大量的工作。根據(jù)文獻(xiàn)[8-10]的結(jié)果, 一般認(rèn)為冗余表示系統(tǒng)對(duì)數(shù)據(jù)進(jìn)行稀疏表示有很好的效果。這些冗余表示系統(tǒng)能夠更好地捕捉不同的圖像特征，并且比正交基具有更佳的稀疏逼近能力。冗余性也使得表示系統(tǒng)的設(shè)計(jì)和訓(xùn)練更加靈活，另外可以帶來(lái)更為魯棒的圖像表示。

本文提出了一種簡(jiǎn)單有效的基于冗余系統(tǒng)的稀疏信號(hào)恢復(fù)模型，求解從部分傅里葉系數(shù)中恢復(fù)原始信號(hào)的問(wèn)題。該方法將L0模與正交投影技術(shù)相結(jié)合，約束優(yōu)化模型的目標(biāo)函數(shù)可以有效利用原始數(shù)據(jù)在給定冗余表示系統(tǒng)中的稀疏性。我們選擇緊小波框架來(lái)對(duì)所恢復(fù)信號(hào)進(jìn)行稀疏逼近，并稱(chēng)該方法為投影迭代硬閾值算法(pIHTA)。L0模的鄰近算子和正交投影技術(shù)被用來(lái)求解算法所產(chǎn)生的子問(wèn)題。由于每個(gè)子問(wèn)題都有閉式解，使得整個(gè)算法簡(jiǎn)潔高效。我們進(jìn)一步證明了pIHTA算法的全局收斂性，并通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證了其有效性。

1 相關(guān)數(shù)學(xué)背景

設(shè)D∈Rm×n(m≥n)，如果?a,b>0，滿(mǎn)足

則稱(chēng)D的行向量是Rn空間的一組框架，當(dāng)a=b時(shí)，稱(chēng)其為緊框架。進(jìn)一步，如果D的列向量規(guī)范正交，即DTD=I，則稱(chēng)D的行向量是Rn空間的一組緊小波框架?？蚣芫哂型陚湫裕覀兛梢杂闷鋵?duì)Rn空間中的任意信號(hào)進(jìn)行表示。由于冗余性的存在，信號(hào)表示方法并不唯一。稀疏表示旨在從所有表示方法中尋找所用非零元系數(shù)最少的一種，即最稀疏的一種表示方法。

本文主要考慮緊小波框架，設(shè)x∈Rn，則α=Dx稱(chēng)為x的小波框架系數(shù)，用其可精確的重建信號(hào)x=DTα。

設(shè)ψ：Rn→(-∞,+∞]是一個(gè)正常的下半連續(xù)函數(shù)，α>0，給定x∈Rn,則由ψ定義的臨近算子proxαψ:Rn→Rn定義為

本文所考慮的不完整傅里葉數(shù)據(jù)的觀測(cè)模型可描述為

y=PFx+η。

(1)

式中：F∈Cn×n表示離散傅里葉變換矩陣；P∈Rd×n(d

2 投影迭代硬閾值算法

由于d

(2)

由于L0模的非凸不連續(xù)性，加之與緊小波框架變換D的復(fù)合，模型(2)的數(shù)值求解異常困難。已有文獻(xiàn)通常將L0模松弛為L(zhǎng)1模，進(jìn)而采用相對(duì)成熟的凸優(yōu)化方法求解[1-3]。但L1模僅是L0模的凸近似，當(dāng)二者之間的等價(jià)性條件不滿(mǎn)足時(shí)，由此導(dǎo)致的模型誤差不可避免。事實(shí)上，采用L1模度量稀疏性通常會(huì)給恢復(fù)信號(hào)帶來(lái)偏差[11]。

本文首先基于約束條件，將L0模和D的復(fù)合運(yùn)算進(jìn)行分離，在緊小波框架變換域進(jìn)行求解，給出模型(2)的如下的等價(jià)模型

(3)

式中：α是x的小波框架系數(shù)；R(D)表示D的值域。顯然模型(2)和(3)的解具有如下關(guān)系：x=DTα。

為了處理模型(3)中的約束條件，我們引入指示函數(shù)lR(α)，

則模型(3)可等價(jià)變形為:

(4)

下面將著重討論等價(jià)模型(4)的數(shù)值求解。為此我們令

顯然g(α)非凸不連續(xù)，f(α)連續(xù)可微且其梯度具有常數(shù)為1的Lipschitz連續(xù)性，即：

模型(4)可由臨近梯度算法求解[6]，其迭代格式為

αk+1=proxγg(αk-γ?f(αk))=

+

(5)

式中γ為下降步長(zhǎng)。雖然公式(5)是一種相對(duì)簡(jiǎn)單的迭代格式，但約束條件α∈R(D)使得αk+1并無(wú)閉式解，導(dǎo)致算法不高效。為了解決該問(wèn)題，類(lèi)似于文獻(xiàn)[12]中的策略，我們采用交替優(yōu)化技術(shù)，首先處理L0模對(duì)應(yīng)項(xiàng)，然后再考慮約束條件，對(duì)應(yīng)的迭代格式為

(6)

式中PR(D)是空間R(D)上的正交投影算子。

下面證明公式(6)中兩個(gè)子問(wèn)題均具有閉式解。對(duì)于第一個(gè)子問(wèn)題由鄰近算子的定義可知

(7)

對(duì)于第二個(gè)子問(wèn)題，由于DTD=I，顯然

注意到

?f(α)=DFTPT(PFDTα-y)，

(8)

該迭代格式具有閉式解，其計(jì)算僅涉及傅里葉變換，緊小波框架變換和硬閾值算子。每步迭代過(guò)程中的信號(hào)可通過(guò)

給出。我們稱(chēng)該格式對(duì)應(yīng)的算法為投影迭代硬閾值算法。

3 收斂性分析

對(duì)于向量α∈Rd,我們用N(α)表示α中零元素的索引集，即

N(α)∶={i|(α)i=0,1≤i≤d}。

我們引入函數(shù)

并且令

易證u(α)連續(xù)可微，并且

?u(α)=DFTPT(PFDTα-y)+

為了估計(jì)?u(α)的Lipschitz常數(shù)，引入常數(shù)

我們有下面的引理

式中λi(A)表示A的第i個(gè)特征值，由于F的正交性和D的列規(guī)范正交性，因此

故

故

綜上可得

(9)

由公式(9)可知，對(duì)?n成立

上式結(jié)合結(jié)論(1)可知結(jié)論(2)成立。

由結(jié)論(2)，?K，對(duì)?k>K，有

則Λ是閉凸集，記PΛ為Λ上的投影算子，當(dāng)k>K時(shí)，我們有

4 數(shù)值實(shí)驗(yàn)

本節(jié)我們用不同的數(shù)據(jù)集，欠采樣模式和緊小波框架來(lái)進(jìn)行數(shù)值實(shí)驗(yàn)，展示pIHTA方法比之前較先進(jìn)的pFISTA方法[12]所恢復(fù)的圖像有更好的品質(zhì)，更快的收斂速度。pFISTA方法在模型(4)中采用L1模度量所恢復(fù)信號(hào)在緊小波框架域的稀疏性，因而是凸優(yōu)化問(wèn)題，便于進(jìn)行理論分析和數(shù)值求解。pFISTA方法是求解相應(yīng)凸優(yōu)化模型的一種快速算法，與本文方法的主要區(qū)別在于迭代過(guò)程中采用軟閾值算子更新系數(shù)α。此外，我們采用相對(duì)誤差(RLNE)：

我們需要設(shè)置幾個(gè)算法參數(shù)。首先根據(jù)多次實(shí)驗(yàn)以及之前成果的經(jīng)驗(yàn)，我們令γ=1，令最大權(quán)重λmax為1，λmin為0.000 1。運(yùn)行算法包括內(nèi)循環(huán)和外循環(huán),內(nèi)循環(huán)更新正則化參數(shù)λ，其收縮率為0.2。我們定義

并基于該相對(duì)誤差設(shè)置內(nèi)、外循環(huán)的停機(jī)準(zhǔn)則。當(dāng)Diff≤10-3時(shí)，λ就會(huì)根據(jù)收縮率變化，當(dāng)Diff≤10-5，算法迭代終止。

平移不變離散小波變換(SIDWT)[13-14]，又稱(chēng)非抽取小波，因?yàn)榭梢暂^好地抑制重建偽影，一般作為典型的緊小波框架用在圖像信號(hào)的恢復(fù)中，我們用具有四個(gè)分解層的Daubechies小波作為SIDWT。

本實(shí)驗(yàn)所用大腦圖像切片(見(jiàn)圖1(a))是3T西門(mén)子掃描儀用加權(quán)渦輪自旋回聲序列掃描一個(gè)健康志愿者得到的。為了增加不同方法對(duì)比結(jié)果的說(shuō)服力，我們同時(shí)選取圖像處理領(lǐng)域中的標(biāo)準(zhǔn)圖像之一“Cameraman”(見(jiàn)圖1(b))進(jìn)行數(shù)值實(shí)驗(yàn)。我們將根據(jù)圖2中的掩模進(jìn)行采樣(白黑像素暗示著采樣和未采樣)。

圖1 測(cè)試圖像

圖2 30%采樣的二維高斯掩模

圖3是在30%采樣率的二維高斯掩模下分別使用pFISTA方法[12]和本文方法恢復(fù)大腦圖片和Cameran圖像的結(jié)果。如果僅從圖像恢復(fù)視覺(jué)效果來(lái)比較，兩種方法差別不大，我們重點(diǎn)關(guān)注兩種算法所恢復(fù)圖像的量化誤差比較以及算法的運(yùn)行時(shí)間。如圖4和表1所示，在相同參數(shù)設(shè)置下，兩組圖像的數(shù)據(jù)實(shí)驗(yàn)均顯示pIHTA方法的重構(gòu)誤差和運(yùn)行時(shí)間(單位：s)優(yōu)于pFISTA方法。

圖3 不同方法的重構(gòu)圖像

圖4 不同方法的相對(duì)誤差和時(shí)間對(duì)比曲線結(jié)果(平移不變離散小波)

表1 兩種方法在不同小波框架下的誤差及運(yùn)行時(shí)間比較

為進(jìn)一步驗(yàn)證本文方法的效果，我們選擇輪廓波(Contourlet)緊小波框架進(jìn)行同樣的數(shù)值實(shí)驗(yàn)[15-16]。在相同的參數(shù)設(shè)置下，運(yùn)行時(shí)間和RLNE的實(shí)驗(yàn)結(jié)果見(jiàn)表1，同樣表明了本文算法的優(yōu)越性。

5 結(jié)語(yǔ)

本文提出了一種投影迭代硬閾值算法來(lái)求解基于不完全傅里葉觀測(cè)數(shù)據(jù)的信號(hào)恢復(fù)問(wèn)題。相應(yīng)的數(shù)學(xué)模型采用L0模度量所恢復(fù)信號(hào)在緊小波框架域的稀疏性。投影迭代硬閾值方法可以有效克服L0稀疏優(yōu)化模型的NP難問(wèn)題，同時(shí)迭代過(guò)程的計(jì)算簡(jiǎn)潔高效。我們證明了該數(shù)值算法所產(chǎn)生序列的全局收斂性。數(shù)值結(jié)果表明與pFISTA方法相比，本文算法具有更好的信號(hào)恢復(fù)效果，收斂速度也相對(duì)更快。本文方法目前尚局限于從不完整的正交變換系數(shù)中恢復(fù)原始信號(hào)，我們將把該方法推廣到一般的線性變換情形，用來(lái)求解諸如壓縮感知、圖像復(fù)原等更廣泛的信號(hào)恢復(fù)問(wèn)題。