一種慣性交替極小化算法及其應(yīng)用

吳雙雙，唐玉超

(南昌大學(xué)數(shù)學(xué)系，江西 南昌 330031)

1 引言與研究背景

在本文中，我們主要考慮以下形式的凸極小化問(wèn)題：

(1)

其中H1和H2是實(shí)Hilbert空間，f:H1→(-∞,+∞]和g:H2→(-∞,+∞]是正則下半連續(xù)凸函數(shù)，G是實(shí)Hilbert空間，A:H1→G和B:H2→G是非零有界線性算子，以及b∈G。我們假設(shè)f是σ-強(qiáng)凸的，其中σ>0。下面，我們介紹幾個(gè)信號(hào)和圖像處理中的具體問(wèn)題，可以看作是問(wèn)題(1)的特例。

例1(全變分圖像去噪問(wèn)題)全變分圖像去噪問(wèn)題最早由Rudin、Osher和Fatemi[1]提出，其形式如下：

(2)

其中u∈Rm×n是觀測(cè)到的噪聲圖像，‖x‖TV表示全變分，λ>0是正則參數(shù)。(2)通常被稱(chēng)為ROF模型。因?yàn)槟Ｐ椭械娜兎帧瑇‖TV可以表示為凸函數(shù)φ和一階差分算子L:Rm×n→Rm×m×Rm×n的組合，即‖x‖TV=φ(Lx)[2-3]。因此，(2)可以表示為下面的形式：

(3)

如果令y=Lx，則(3) 可以寫(xiě)成如下形式：

s.t.Lx-y=0

(4)

例2(具有約束的全變分圖像去噪問(wèn)題) 考慮到數(shù)字圖像中像素值往往分布在[0,1]或[0,255]之間。設(shè)C是一個(gè)非空閉凸集，考慮如下具有約束的ROF模型[4-6]：

s.t.x∈C

(5)

通過(guò)使用示性函數(shù)δC(x)，以及‖x‖TV=φ(Lx)， (5)等價(jià)于：

(6)

s.t.Lx-y=0

(7)

這也是(1)的一個(gè)特例。除此之外，還有許多問(wèn)題是(1)的特例，例如閉凸集交集上的投影問(wèn)題[7-8]，強(qiáng)凸魯棒主成分分析(RPCA)問(wèn)題[9-10]和資源分配問(wèn)題[11]等。

交替方向乘子法(ADMM)是求解(1)的一種常用方法，最早可以追溯到Gabay和Mercier[12]、Glowinski和Marroco[13]以及Gabay[14]等人的工作。眾所周知，ADMM等價(jià)于Douglas-Rachford分裂算法[15]應(yīng)用于(1)的對(duì)偶問(wèn)題。為了充分利用(1)中f(x)的強(qiáng)凸性，基于鄰近梯度算法(也稱(chēng)為向前向后分裂算法)，Tseng[16]提出了一種交替極小化算法(AMA)求解問(wèn)題(1)的對(duì)偶問(wèn)題，具體算法如下：

(8)

為了加速AMA(8)，Goldstein等[17]基于快速迭代收縮閾值算法 (FISTA)[18]提出了一種求解(1)的對(duì)偶問(wèn)題的算法，即Fast AMA：

(9)

由于Fast AMA (9) 是由FISTA推導(dǎo)建立的，其步長(zhǎng)參數(shù)范圍比原始AMA (8)小，而文獻(xiàn)[17]并沒(méi)有證明序列{xk}和{yk}的收斂性。因此，本文的主要目的是提出一種用于求解(1)的慣性交替極小化算法，該算法源于慣性鄰近梯度算法。我們不僅分析所提出算法的收斂性，而且證明所提出的算法收斂到原始問(wèn)題(1)和對(duì)偶問(wèn)題的最優(yōu)解。作為應(yīng)用，我們將所得算法用于求解一類(lèi)復(fù)合凸極小化問(wèn)題。為驗(yàn)證所提出算法的有效性和優(yōu)越性，我們考慮具有約束的全變分圖像去噪模型，該模型包含非空閉凸集約束和全變分正則項(xiàng)。

本文的主要貢獻(xiàn)：(ⅰ)提出一種慣性的交替極小化算法求解凸極小化問(wèn)題(1)，其中AMA(8)可以看作本文所提算法的特例。(ⅱ)證明所提出算法同時(shí)收斂到問(wèn)題(1)和其對(duì)偶問(wèn)題的最優(yōu)解。(ⅲ)應(yīng)用所提算法求解具有約束的全變分圖像去噪模型，數(shù)值結(jié)果表明了所提算法的優(yōu)越性。

本文具體安排如下：在第2節(jié)，我們回顧凸分析中的一些基本的定義和重要結(jié)論。在第3節(jié)，我們介紹所提算法并證明其收斂性。在第4節(jié)，我們利用提出的算法來(lái)解決一類(lèi)復(fù)合凸優(yōu)化問(wèn)題。在第 5節(jié)，通過(guò)數(shù)值實(shí)驗(yàn)說(shuō)明所提出算法的性能。最后，我們給出全文小結(jié)。

2 預(yù)備知識(shí)

C上的投影用PC表示，定義為PC(x)=argminy∈C‖x-y‖。

定義1[20]設(shè)f:H→(-∞,+∞]，

(ⅰ)如果

f(λx+(1-λ)y)≤λf(x)+(1-λ)f(y),?x,y∈domf,λ∈[0,1]

那么f是凸的。

f在x處的次微分定義為:

?f(x)≡{v∈H|f(y)≥f(x)+〈v,y-x〉,?y∈H}

注意f是凸的就意味著?f是單調(diào)的，即

〈x-y,u-v〉≥0,?x,y∈H,?u∈?f(x),?v∈?f(y)

此外，如果f是σ-強(qiáng)凸的，則?f是σ-強(qiáng)單調(diào)的，即?f-σI是單調(diào)的。

定義2[20]設(shè)f∈Γ0(H)，則f的Fenchel共軛定義為

由f∈Γ0(H)，可以得到(?f)-1=?f*。

定義3[20](鄰近算子)設(shè)λ>0，f∈Γ0(H)，λf的鄰近算子定義為：

如果f(x)=δC(x)，那么proxλδC(x)=PC(x)。

Moreau[21]首次提出了鄰近算子的概念，它在設(shè)計(jì)和解決凸優(yōu)化問(wèn)題的鄰近算法中起著非常重要的作用。下面的引理展示了如何利用一個(gè)正則下半連續(xù)凸函數(shù)f的凸共軛f*或其逆來(lái)計(jì)算它本身的鄰近算子。

引理1[20](Moreau's decomposition)設(shè)λ>0，f∈Γ0(H)，則

定義4[20](Moreau envelope)設(shè)f∈Γ0(H)，λ>0，則Moreau envelope可以表示為：

下面的引理表明Moreau envelope在H上是可微的。

2.1 慣性鄰近梯度算法

考慮下面凸極小化問(wèn)題：

(10)

其中f:H→R是具有L-Lipschitz連續(xù)梯度的凸可微函數(shù)，L>0，g:H→(-∞,+∞]是一個(gè)正則的、下半連續(xù)凸函數(shù)。鄰近梯度算法 (PGA) 是求解(10)的一種常用的方法。PGA 可以看作是在凸極小化問(wèn)題中向前-向后分裂算法[22-24]的一種特殊情況。此外，慣性鄰近梯度算法是一種有效的加速鄰近梯度算法被廣泛應(yīng)用于求解問(wèn)題(10)。下面，我們給出關(guān)于慣性鄰近梯度算法的收斂性定理，具體證明可以參考文獻(xiàn)[25]。

定理1設(shè)g∈Γ0(H)，f:H→R是具有L-Lipschitz連續(xù)梯度的凸可微函數(shù)，其中L>0。假設(shè) argmin(f+g)≠?，x0∈H，定義

(11)

如果{αk}和{γk}滿足下列條件之一：

則以下結(jié)論成立：

(b){xk}弱收斂到argmin(f+g)的一個(gè)極小點(diǎn)。

2.2 最優(yōu)性條件

對(duì)于前面我們所考慮的凸極小化問(wèn)題(1)，它的拉格朗日函數(shù)定義為

L(x,y,p)=f(x)+g(y)+〈p,b-Ax-By〉

(12)

我們通過(guò)極小化關(guān)于x和y的拉格朗日函數(shù)(12)來(lái)獲得對(duì)偶問(wèn)題，如下所示：

(13)

其中f*和g*分別是f和g的Fenchel共軛。如果

L(x*,y*,p)≤L(x*,y*,p*)≤L(x,y,p*),?(x,y,p)∈H×G×K

則稱(chēng)(x*,y*,p*)是拉格朗日函數(shù)L的鞍點(diǎn)。眾所周知，當(dāng)且僅當(dāng) (x*,y*)是(1)的最優(yōu)解，p*是(13)的最優(yōu)解時(shí)，(x*,y*,p*)是拉格朗日函數(shù)L的鞍點(diǎn)，并且原始問(wèn)題(1) 的最優(yōu)值和對(duì)偶問(wèn)題(13)的最優(yōu)值是相同的。(x*,y*,p*) 通常稱(chēng)為(1)和(13) 的原始-對(duì)偶最優(yōu)解。

引理3[11](x*,y*,p*)是拉格朗日函數(shù)L的鞍點(diǎn)，當(dāng)且僅當(dāng)

A*p*∈?f(x*)

B*p*∈?g(y*)

A*x*+B*y*=b

(14)

最后，在本文中我們將充分利用以下余弦等式：

2〈a-b,c-a〉=‖b-c‖2-‖a-b‖2-‖a-c‖2,?a,b,c∈H

(15)

3 慣性交替極小化算法收斂性分析

在本節(jié)中，我們提出一種用于求解(1)的慣性交替極小化算法，它可以看作是原始AMA算法(8)的推廣。下面我們給出慣性交替極小化算法的具體格式。

算法1慣性的交替極小化算法

當(dāng)αk=0時(shí),所提出的IAMA退化成原始的AMA(8)。下面，我們證明所提的慣性交替極小化算法的收斂性。

定理2考慮問(wèn)題(1)，假設(shè)拉格朗日函數(shù)(12)的鞍點(diǎn)集是非空的，{xk},{yk}和{pk} 是由算法1生成的序列，如果{αk}和{γk}滿足下列條件之一：

(16)

(17)

則存在拉格朗日函數(shù)L的一組鞍點(diǎn)(x*,y*,p*)，使得：

(a)pk?p*；

(b)xk→x*；

(c)Byk→By*；此外，如果存在β>0，使得B*BβI，則yk→y*；

證明(a)由引理4可知，f*(A*p) 是連續(xù)可微的，并且具有Lipschitz連續(xù)梯度。在對(duì)偶問(wèn)題(13)中，極小化的是一個(gè)光滑函數(shù)和一個(gè)凸函數(shù)的和，因此我們可以應(yīng)用慣性鄰近梯度算法來(lái)求解(13)，即：

(18)

(19)

即

(20)

將h(p)代入(18)的第二個(gè)式子，有

(21)

由(21)可得到其一階優(yōu)化條件為

(22)

設(shè)yk+1∈?g*(B*pk+1),我們有

0∈?g(yk+1)-B*pk+1

(23)

以及

(24)

這是我們所提出算法的最后一步。由(23)和(24)，我們可以得到

(25)

這是慣性交替極小化算法的第三步。根據(jù)定理 1，我們可以得到pk?p*。

(b),(c)由一階優(yōu)化條件和算法可得到

(26)

和

(27)

根據(jù)?f強(qiáng)單調(diào)和?g單調(diào)，則有

(28)

和

(29)

將(28)和(29)相加，可以得到

(30)

根據(jù)pk+1的定義，以及Ax*+By*=b，可以得到

(31)

將(31)代入(30)，可以得到

(32)

對(duì)(32)利用余弦等式，我們有

(33)

接著在上述不等式兩邊乘以2γk，經(jīng)過(guò)移項(xiàng)和整理可以得到

(34)

(35)

設(shè)φk=‖pk-p*‖2，上式可以重新改寫(xiě)為

(36)

將(36)重新整理，可以得到

(37)

將(34)再次進(jìn)行移項(xiàng)，可以得到

(38)

對(duì)(38)兩邊同時(shí)進(jìn)行求和，則

(39)

在(39)兩邊取極限，即當(dāng)k→+∞時(shí)，我們能夠得到

即

xk→x*,Byk→By*=b-Ax*

如果β>0，使得B*BβI，那么yk→y*。

(d)下面我們證明目標(biāo)函數(shù)值f(xk+1)+g(yk+1)收斂于(1)的最優(yōu)解。由(12)所定義的拉格朗日函數(shù)L我們可以得到：

L(xk+1,yk+1,p*)≥L(x*,y*,p*)

(40)

把它展開(kāi)來(lái)，即：

f(xk+1)+g(yk+1)≥f(x*)+g(y*)+〈p*,Axk+1+Byk+1-b〉

(41)

(42)

由算法1的一階優(yōu)化條件我們可以得到

(43)

由次微分的定義，可以得到：

(44)

把上面兩個(gè)式子相加得到f(x*)+g(y*)≥f(xk+1)+g(yk+1)+〈A*(pk+αk(pk-pk-1)),x*-xk+1〉+〈pk+1,By*-Byk+1〉

(45)

(46)

最后，結(jié)合(42)和(45)，可以得出

(47)

注1在定理 2中，我們證明了慣性交替極小化算法在兩種條件下的收斂性。

(ⅱ)在第二種情況下，為了方便，我們固定步長(zhǎng)參數(shù)γk，記為γ，根據(jù)文獻(xiàn)[25]中的定理 2， (17)可以寫(xiě)為

因?yàn)樾蛄衶αk}是非遞減的，所以0≤αk≤α(γ)，其中

α(γ)=1+

(48)

在本文中，我們?nèi)ˇ?10-6。

4 應(yīng)用

在本節(jié)中，我們應(yīng)用所提出的算法1來(lái)求解一類(lèi)復(fù)合凸極小化問(wèn)題，此類(lèi)問(wèn)題在信號(hào)和圖像處理中具有廣泛的應(yīng)用。具體來(lái)說(shuō)，我們考慮如下形式的凸極小化問(wèn)題。

問(wèn)題4.1假設(shè)g∈Γ0(H)、h∈Γ0(G)、b∈H，并且L:H→G是一個(gè)非零有界線性算子?？紤]如下凸極小化問(wèn)題：

(49)

我們對(duì)所考慮的凸極小化問(wèn)題 4.1做以下假設(shè)：(i)，最優(yōu)解存在；(ii)，滿足一定的正則性條件，例如0∈sqri(L(domg)-domh)。 容易看出問(wèn)題 4.1 是問(wèn)題(1) 的一個(gè)特例。下面，我們給出具體的收斂性定理。

定理3在問(wèn)題4.1的設(shè)定中，設(shè)p0∈G，定義

(50)

其中{γk}和{αk}滿足定理2的條件。則下列說(shuō)法成立:

(ⅰ){pk}弱收斂到(51)的對(duì)偶問(wèn)題的一個(gè)最優(yōu)解p*，其對(duì)偶問(wèn)題為：

(51)

(ⅱ)x*=proxg(u+L*p*)是(51)的唯一最優(yōu)解，并且{xk}強(qiáng)收斂到x*。

s.t.Lx-y=0

(52)

易見(jiàn)(52)是(1)的一種特殊情況。因此，我們可以應(yīng)用算法1來(lái)求解(52)。從而我們得到如下迭代格式的算法：

(53)

由xk+1的一階優(yōu)化條件，可以得到

(54)

再由(53)中yk+1的定義和引理1，可以得到

yk+1=

(55)

將(55)代入到(53)的pk+1中，可以得到

(56)

即有(50)。

(ⅰ)根據(jù)定理2，我們知道{xk}強(qiáng)收斂到x*,{(yk,pk)} 弱收斂到 (y*,p*)，其中(x*,y*,p*) 是下面拉格朗日函數(shù)的鞍點(diǎn)：

L(x,y,p)=f(x)+h(y)+〈p,-Lx+y〉

(57)

由引理3，我們得到

(58)

根據(jù)(13)我們可以得到(52)的對(duì)偶問(wèn)題是：

(59)

(ⅱ)由(58) 知，x*=proxg(u+L*p*)是(49)的唯一解，且{xk}強(qiáng)收斂到x*。

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

在本節(jié)中，我們應(yīng)用所提算法求解具有約束的全變分圖像去噪模型(5)，并與其他算法進(jìn)行比較。所有實(shí)驗(yàn)均在Windows 7和MATLAB R2016a下進(jìn)行，該筆記本電腦的配置為 Intel Core i7-6700 CPU @3.40GHZ和8G RAM。在下文中，我們將具有約束的全變分圖像去噪模型(5)稱(chēng)為C-L2-TV。通常存在兩種類(lèi)型的全變分：一種是各向同性全變分 (ITV)，另一種是各向異性全變分 (ATV)。它們都可以用凸函數(shù)φ和一階差分算子L的組合來(lái)表示，即‖x‖TV=φ(Lx)。本文所提出的算法可應(yīng)用于 ITV和 ATV。在下面的實(shí)驗(yàn)中，我們只使用ATV作為(5)中的正則函數(shù)，類(lèi)似地可以處理ITV的情況。我們定義約束集合C={x|0≤x≤255}。

我們用峰值信噪比 (PSNR) 和結(jié)構(gòu)相似性指數(shù) (SSIM) [26]來(lái)衡量恢復(fù)圖像的質(zhì)量，其定義如下：

和

實(shí)驗(yàn)測(cè)試圖像包括“Cameraman” 、 “Building”、“Boat”和“Goldhill”，如圖1所示。

(a) 256×256 “Cameraman”圖像(b) 493×517 “Building”圖像 (c) 512×512 “Boat”圖像 (d) 512×512 “Goldhill”圖像

表1 不同圖像在不同噪聲水平下的最優(yōu)正則參數(shù)選擇Table 1 Optimal regularization parameter selection for different images under different noise levels

我們同時(shí)考慮了兩種不同慣性參數(shù)條件下的IAMA，分別稱(chēng)為IAMA_1和IAMA_2。對(duì)原始圖像分別添加均值為0和標(biāo)準(zhǔn)方差分別為σ=15,25和50的高斯噪聲。為了獲得最佳的去噪圖像，我們調(diào)整了正則化參數(shù)，所得數(shù)值結(jié)果見(jiàn)表2。

在第一個(gè)實(shí)驗(yàn)中，我們比較本文提出的慣性AMA算法與原始的AMA算法在不同步長(zhǎng)參數(shù)下的恢復(fù)效果。對(duì)于模型(5)，σ=1，‖L‖2=8，由注1，我們有

為了直觀比較，我們給出α(γ)的圖像，見(jiàn)圖2。

圖2 由γ決定的慣性參數(shù)αk的上界Figure 2 The upper bound of the inertial parameter αk determined by γ

接下來(lái)我們?cè)O(shè)置了四種不同的步長(zhǎng)參數(shù)來(lái)進(jìn)行試驗(yàn)，γk分別選取：γk=0.1,γk=0.15,γk=0.2,γk=0.225。在IAMA_2中，根據(jù)注1我們?nèi)。害羕=0.26,αk=0.207,αk=0.132,αk=0.078。這里的測(cè)試圖像為“Building”，我們對(duì)原始圖像添加均值為0 和標(biāo)準(zhǔn)方差為σ=50的高斯噪聲。

所得數(shù)值結(jié)果如表2所示，由表中數(shù)據(jù)我們可以看出來(lái)，在不同的步長(zhǎng)參數(shù)下，iAMA_1和原始的AMA在峰值信噪比 (PSNR)、結(jié)構(gòu)相似性指數(shù) (SSIM) 以及迭代步數(shù)的數(shù)值結(jié)果都是一致的。在給定停止條件下，隨著γk在取值的增大，3種算法所需迭代次數(shù)逐漸減少。同時(shí)，IAMA_2算法比原始的AMA算法收斂稍快。

表2 慣性AMA和原始AMA在不同參數(shù)關(guān)于PSNR、SSIM 和 Iter的數(shù)值結(jié)果Table 2 Numerical results of the proposed inertial AMA and the AMA in terms of the PSNR,SSIM and Iter with different parameters

在第2個(gè)實(shí)驗(yàn)中，我們測(cè)試了慣性AMA和原始AMA在對(duì)于不同圖像且在不同噪聲下的恢復(fù)效果。我們統(tǒng)一設(shè)置γk=0.1。在IAMA_1中，我們按照注1對(duì)慣性參數(shù)來(lái)進(jìn)行更新。在IAMA_2中，我們?cè)O(shè)定慣性參數(shù)αk=0.26。

表3展示了4幅圖像在不同噪聲水平下去噪效果的數(shù)值結(jié)果。由表中數(shù)據(jù)可以看出來(lái)，IAMA_1和原始的AMA所選取的圖像在不同的噪聲水平下的峰值信噪比(PSNR)、 結(jié)構(gòu)相似性指數(shù)(SSIM)以及迭代步數(shù) (Iter) 幾乎都是一樣的。對(duì)于所有的測(cè)試圖像，在不同噪聲水平下，IAMA_2的迭代步數(shù)都是要比原始的AMA迭代步數(shù)少。在4幅測(cè)試圖像中，噪聲水平比較大的時(shí)候，即σ=50時(shí)，IAMA_2的恢復(fù)效果都要比原始的AMA恢復(fù)效果稍好。在“Cameraman”圖像中，對(duì)于設(shè)定的3種噪聲水平，IAMA_2 的恢復(fù)效果都要比原始的AMA恢復(fù)效果稍好?？傮w來(lái)說(shuō)，我們所提出的算法比原始的AMA算法具有一定的優(yōu)越性。接下來(lái)，我們展示“Cameraman”圖像在不同噪聲水平下，本文所提算法與原始AMA算法的恢復(fù)效果，如圖3所示。

表3 所提算法在兩種不同的慣性參數(shù)設(shè)置下與AMA的數(shù)值結(jié)果Table 3 Numerical results of the AMA and the proposed algorithm under two different inertial parameter

(a)σ=15 (b)σ=25 (c)σ=50

(d)AMA (e)AMA (f)AMA

(g)IAMA_1 (h)IAMA_1 (i)IAMA_1

(j)IAMA_2 (k)IAMA_2 (l)IAMA_2

6 結(jié)論

為解決具有線性等式約束的兩塊可分離凸極小化問(wèn)題(1)，我們提出了一種慣性交替極小化算法，所提出的算法推廣了 Tseng[16]提出的經(jīng)典交替極小化算法。在實(shí)Hilbert空間中，我們證明了所提出的算法在兩種不同慣性參數(shù)條件下的收斂性。進(jìn)一步，我們應(yīng)用所提出的算法求解一類(lèi)復(fù)合凸極小化問(wèn)題，該問(wèn)題在圖像處理，特別是圖像去噪中具有廣泛的應(yīng)用。為了驗(yàn)證所提出算法的有效性和優(yōu)越性，我們應(yīng)用所提出的算法求解具有約束的全變分圖像去噪模型。數(shù)值實(shí)驗(yàn)表明慣性交替極小化算法比原始交替極小化算法的迭代步數(shù)少。