吳慶豐
(淮北師范大學 數(shù)學科學學院,安徽 淮北 235000)
不等式約束優(yōu)化問題是約束優(yōu)化問題中一類重要問題,一種經(jīng)典的解法是采用SUMT(sequential un?constrained minimization technique)技術(shù),即把約束條件轉(zhuǎn)化為某種懲罰函數(shù)加到目標函數(shù)中去,從而將約束優(yōu)化轉(zhuǎn)化為一系列無約束優(yōu)化問題的求解.本文將采用一種新的內(nèi)點算法求解,并與傳統(tǒng)內(nèi)點障礙函數(shù)法進行比較.下面介紹求解凸優(yōu)化問題的IPA(inerior point algorithm)方法.
文獻[1]和文獻[2]中分別給出Bregman函數(shù)和Bregman-Legendre函數(shù)的定義如下:
設(shè)S為開凸集,S?Rn,對可微函數(shù)定義Bregman距離Dg(x,y)=g(x)-g(y)-<?g(y),x-y>.稱g為區(qū)域S上的Bregman函數(shù),如果下列條件被滿足:
(B1)g于連續(xù)且嚴格凸;
(B2)g于S連續(xù)可微;
(B3)Dg(x,·)及Dg(·,y)的水平集有界,
作為對Bregman 函數(shù)的擴展,在文獻[2]中Bauschke 等引入了Bregman-Legendre 函數(shù)的概念.設(shè)f:RJ→(-∞,+∞)為適當?shù)拈]凸函數(shù),其本質(zhì)區(qū)域定義為這里,我們稱凸函數(shù)為適當?shù)?,如果不存在x使得f(x)=-∞,且存在x滿足f(x)<+∞;稱適當?shù)耐购瘮?shù)f為閉的,如果f是下半連續(xù)的.記Df(·,·):D×intD→[0,+∞),Df(x,z)=f(x)-f(z)-<?f(z),x-z>,這里 intD表示D的內(nèi)部,稱Df(·,·)為Bregman距離.下面給出求解凸優(yōu)化問題的IPA方法.
IPA算法[3]給定x0∈intD,對k=0,1,…,求xk+1滿足,
這里假定每步迭代中求得的xk+1∈intD.
關(guān)于該算法的收斂性,Byrne給出了如下定理.
定理1[3]設(shè)存在.如是被唯一確定的,則由如上IPA算法產(chǎn)生的點列{k} 收斂于;如不能被唯一確定,但可在D中選取,則序列{Dh(,k)}和{DF(,xk)}單減.如果進一步,DF水平集有界,則{xk} 有界,且其任一聚點x*∈滿足f(x*)=f(x).如h或F為 Bregman-Legendre 函數(shù),本質(zhì)區(qū)域為D,則x*∈D,且點列收斂于x*.
注1:這里DF水平集有界,是指對x∈D及 ?a>0,集合{z|DF(,z)≤a} 有界.
下面討論如下不等式約束的凸優(yōu)化問題的求解:
這里f(x)是Rn中連續(xù)可微的凸函數(shù),gi(x),i=1,2,…,m是Rn中連續(xù)可微的凹函數(shù).
對于上述問題,一種傳統(tǒng)的方法是采用對數(shù)障礙函數(shù)法求解,即通過求解一系列無約束凸優(yōu)化問題
來得到(2.1)的解,這里r>0 是罰參數(shù).
記D={x|gi(x)>0,i=1,2,…,m},易知D為開凸集.由于gi(x),i=1,2,…,m是連續(xù)可微的凹函數(shù),可知為D上連續(xù)可微的凸函數(shù),又由于f(x)為連續(xù)可微凸函數(shù),故對r>0,若問題(2.2)存在極小點,則在該點處滿足
即對數(shù)障礙函數(shù)法的每步迭代中,xk+1滿足
在對數(shù)障礙函數(shù)法的計算中,要求罰參數(shù)r單減并趨于0,所以當r充分小時,即使嚴格凸,也可能出現(xiàn)病態(tài)問題,使得求解難度增大.那么可否避免這一問題呢?
假設(shè)Slater條件成立,即D≠?.由于我們只關(guān)心f(x)在D上的取值,不妨假定f(x)在D外的值均為+∞,則f(x)為適當?shù)拈]凸函數(shù).令F(x)=f(x)+rh(x),取h(x)使之滿足IPA 所要求的條件,基于IPA 可以得到如下算法:
求xk+1∈D使之滿足 ?F(xk+1)=r?h(xk)或 ?f(xk+1)+r?h(xk+1)=r?h(xk).
令
容易驗證h(x)滿足IPA所要求的條件.于是,基于IPA可以得到如下算法:
算法1取x0∈D及r>0,對k=0,1,2,…,計算xk+1∈D滿足:
假定每步迭代中滿足條件的xk+1存在.
算法1 需要假定每步迭代中xk+1可以在D中求得,也就是當求得的迭代點xk+1?D時,算法無法實現(xiàn).下面的關(guān)鍵問題是如何求解(2.5)式,并且保證xk+1∈D.
f(x)和h(x)都是x的凸函數(shù),則F(x)也是x的凸函數(shù).設(shè)F(x)是二次可微實函數(shù),給定初始xk∈Rn,Hesse矩陣?2F(xk)正定.我們在xk附近用二次Taylor展開,
取xk+1=x*,令 ?F(x*)=0,得xk+1=xk-[?2F(xk)]-1?F(xk),此即牛頓法迭代公式.由 IPA 算法可知,?F(x*)=r?h(xk),則
由F(x)=f(x)+rh(x)得
此公式與牛頓法公式相似.我們將xk+dk看作(1.1)的一個近似解,這里
令xk+1=xk+dk,如果xk+1∈D,則可以作為下次迭代的初始點.若xk+1?D,則添加步長αk(0<αk<1),使xk+1=xk+αkdk∈D.顯然,當xk+1=xk+dk∈D時直接作為新的迭代點,實際上是默認步長αk=1.
借鑒阻尼牛頓法,加入線搜索技術(shù),搜索準則采用Armijo-Goldstein 準則或Wolfe-Powell準則求步長αk,并且使目標函數(shù)值有一定的下降,當xk+1=xk+αkdk∈D時,步長保持不變;當xk+1=xk+αkdk?D時,減小步長αk直到xk+1=xk+αkdk∈D.在計算搜索方向dk時,增廣目標函數(shù)的Hesse矩陣仍然有可能奇異或者病態(tài).這時,考慮搜索方向結(jié)合最速下降法,即當?2F(xk)不正定時,令Hesse矩陣為單位矩陣.這時,下降方向為采用最速下降方向.下面給出如下算法.
算法2:
步一 給定終止誤差 0≤ε?1,r> 0,0<c<1,0<η<1,初始點x0∈D,k=1.
步二 若‖?f(xk)‖≤ε,則x*=xk,算法終止;否則計算搜索方向,記
步三令xk+1=xk+αkdk,若xk+1∈D,k=k+1,轉(zhuǎn)步二;否則縮小步長,令αk=cαk,直 到xk+1=xk+αkdk∈D,k=k+1,轉(zhuǎn)步二.
下降方向dk結(jié)合最速下降法時可采用原目標函數(shù)負梯度或者增廣目標函數(shù)的負梯度方向,線搜索可以從原目標函數(shù)或者增廣目標函數(shù)考慮.其實計算的最終目的就是要原目標函數(shù)下降,當然又要保證迭代點在可行域內(nèi).顯然,當?shù)c靠近或超出可行域邊界時,增廣目標函數(shù)將趨于無窮大,根據(jù)這一點,可通過適當減少步長控制每一步迭代過程中迭代點仍為內(nèi)點.
上述算法中只需事先給定一個r,不必趨于零,所以避免了傳統(tǒng)內(nèi)點障礙函數(shù)法由于rk→0,(k→∞)時,增廣目標函數(shù)病態(tài)問題的出現(xiàn).線搜索可以從原目標函數(shù)或者增廣目標函數(shù)考慮.從數(shù)值實驗中計算效果來看,線搜索直接考慮使原目標函數(shù)有充分的下降量計算時間更短.其實計算的最終目的就是要原目標函數(shù)下降,當然又要保證迭代點在可行域內(nèi).這里的內(nèi)點算法給出了一種新的搜索方向,計算過程中不需要改變r,加入步長控制技術(shù),使迭代點在可行域內(nèi)部.
這里給出一個不等式約束的凸優(yōu)化的例子,采用算法2進行計算.本例的整體極小點容易看出,這樣是為了方便與用算法2所計算的結(jié)果進行比較.
求解下面優(yōu)化問題.
容易看出此問題的整體極小點為x*=(1,0)T.取初始點x0=(2,1)T,這里采用算法2,利用Matlab編程求解,求解結(jié)果如表1所示.
表1 算法2求解的數(shù)值結(jié)果
r的取值不必很小,從數(shù)值實驗結(jié)果看,計算效果還不錯.實際上,r如果取得太小,計算中就會導致病態(tài)問題的出現(xiàn),這一點正是傳統(tǒng)內(nèi)點障礙函數(shù)法的一個主要缺點.從數(shù)值實驗中計算效果來看,線搜索直接考慮使原目標函數(shù)有充分的下降量計算時間更短.
本文提出的內(nèi)點算法,給出了一種非精確求解IPA的方法,從而便于將算法在計算機上實現(xiàn),并且還解決了IPA對迭代點在可行域內(nèi)部的假定.算法中只需事先給定一個r,不必趨于零,所以避免了傳統(tǒng)內(nèi)點障礙函數(shù)法由于rk→0,(k→∞)時,增廣目標函數(shù)病態(tài)問題的出現(xiàn).這里的內(nèi)點算法給出了一種新的搜索方向,計算過程中不需要改變r,加入步長控制技術(shù),使迭代點在可行域內(nèi)部.關(guān)于算法的收斂性證明將是另一重要內(nèi)容.
[1]BREGMAN L M.The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[J].USSR Computational Mathematics and Mathematical Physics,1967,7(3):200-217.
[2]BAUSCHKE H H,BORWEIN J M.Legendre functions and the method of random Bregman projections[J].Journal of Con?vex Analysis,1997,4(1):27-67.
[3]BYRNE C.Block-iterative interior point optimization methods for image reconstruction from limited data[J].Inverse Prob?lems,2000,16(5):1405-1419.