基于光滑逼近函數的高階牛頓法求解凸二次規劃

雍龍泉, 賈 偉, 黎延海

(1.陜西理工大學數學與計算機科學學院, 漢中 723001； 2.陜西理工大學, 陜西省工業自動化重點實驗室, 漢中 723001)

凸二次規劃(convex quadratic programming,CQP)是非線性規劃中的一類特殊數學規劃問題,在很多方面都有應用,如投資組合、支持向量機等[1-5]。鑒于求解凸二次規劃是模式識別、機器學習等領域的關鍵環節,因此(大規模稀疏)凸二次規劃一直是優化理論的研究熱點。求解二次規劃常見方法有:拉格朗日方法、Lemke方法、內點法、有效集法和橢球算法等[6-11]。自從1984年印度著名數學家Karmarker提出了內點算法以來,由于該算法具有多項式復雜性,而且實際計算效果較為良好,經過多年的研究,凸二次規劃的原對偶內點算法以及核函數內點算法已經成功地被推廣到求解半定優化、互補問題、錐規劃等優化問題。凸二次規劃已經成為運籌學、經濟數學、人工智能、機器學習等領域的基本方法。

盡管內點算法具有多項式迭代復雜性,但是每次計算下降方向時采用牛頓方向:即計算一個齊次線性方程組的非零解。這其中需要計算雅克比矩陣及其逆矩陣,雅克比矩陣的非奇異性保證了齊次線性方程組解(下降方向)的存在唯一性。 而實際問題的雅克比矩陣往往是一個高維的稀疏矩陣,其行列式或特征值難以計算,因此雅克比矩陣的非奇異性難以證明,相應的齊次線性方程組求解也不是太容易,這也是目前內點算法的瓶頸。

為克服各類內點算法中計算瓶頸,即不考慮牛頓方向。給出求解凸二次規劃的新方法：通過將凸二次規劃轉化為非線性方程組,進而采用逼近程度較好(且不易發生溢出的)光滑逼近函數進行處理,得到光滑非線性方程組,再采用三步迭代牛頓法進行求解。

1 光滑逼近函數

定義1給定非光滑函數f(t):R→R,稱光滑函數fμ(t),μ>0為f(t)的光滑逼近函數,如果對任意t∈R,存在κ>0,使得|fμ(t)-f(t)|≤κμ,?μ>0，如果κ不依賴于t,則稱fμ(t)為f(t)的一致光滑逼近函數。

一致光滑逼近函數可以分為從上方一致逼近和從下方一致逼近[12]。文獻[13]給出絕對值函數φ(t)=|t|的光滑函數:

(1)

(2)

(3)

3個一致光滑逼近函數的性質如下。

(4)

(3)對于?t∈R,有

(5)

圖1 μ=0.4時與φ(t)=|t|的圖像Fig.1 Graph of and φ(t)=|t| with μ=0.4

圖2 μ=0.2時與φ(t)=|t|的圖像Fig.2 Graph of and φ(t)=|t| with μ=0.2

2 凸二次規劃光滑處理

考慮凸二次規劃及對偶二次規劃(dual quadratic programming,DQP)。

(1)CQP：

(6)

(2)DQP:

(7)

式中：Q∈Rn×n為對稱的半正定矩陣；c∈Rn；z∈Rn;A∈Rm×n;b∈Rm;w∈Rn;y∈Rm;m≤n,且rank(A)=m,其最優性條件為

(8)

令z=|x|-x,w=|x|+x,則最優性條件等價于非線性方程組:

(9)

記向量絕對值函數φ(x)=|x|=(|x1|,|x2|,…,|xn|)T的每一個分量為φ(xi)=|xi|,i=1,2,…,n。

對每一個φ(xi),采用定理1中的光滑函數

(10)

進行光滑化處理,這樣就得到向量絕對值函數φ(x)的光滑逼近函數為

φμ(x)=[φμ(x1),φμ(x2),…,φμ(xn)]T

(11)

且每一個分量φμ(xi)滿足0<φμ(xi)-φ(xi)≤μ,i=1,2,…,n。

從而有

(12)

通過光滑處理,問題[式(9)]轉化為光滑方程組:

0

(13)

定理2對任意的μ>0,Fμ(x,y)一致光滑逼近F(x,y)。

證明任意的μ>0,

Fμ(x,y)-F(x,y)=

(14)

引理1 設矩陣A∈Rm×n,rank(A)=m,D=diag(d1,d2,…,dn), |di|<1,i=1,2,…,n。Q∈Rn×n為對稱半正定矩陣,I表示單位矩陣,則矩陣

(15)

是非奇異的。

證明令X=-(I+D)(I-D)-1-Q,其中，-X為正定矩陣,故X可逆。則有

(16)

故G的非奇異性等價于-AX-1AT的非奇異性。

由于-X=(I+D)(I-D)-1+Q為正定矩陣,結合rank(A)=m可知，A(-X-1)AT為正定矩陣,從而-AX-1AT非奇異,也即G非奇異。

記Fμ(x,y)的Jacobian矩陣F′μ(x,y), 簡單運算可得

(17)

式(17)中：

(18)

定理3設Q∈Rn×n為對稱半正定矩陣,A∈Rm×n,rank(A)=m,則對?μ>0,Fμ(x,y)的Jacobian矩陣F′μ(x,y)非奇異。

證明根據定理1,則有

(19)

結合引理1可知,F′μ(x,y)非奇異。凸二次規劃的最優性條件[式(18)]已經轉化為光滑非線性方程組[式(13)],給出求解光滑非線性方程組[式(13)]的高階牛頓法。

3 高階牛頓法

(1)輸入。矩陣Q∈Rn×n,A∈Rm×n, rank(A)=m,b∈Rm,c∈Rn；設定常數μ和終止條件ε；取x(0)∈Rn,y(0)∈Rm；令p=0。

(2)開始計算。

X(p)=(x(p),y(p))

(20)

(21)

Fμ(X(p))=

(22)

(23)

F′μ(X(p))=

(24)

(25)

V(p)=X(p)-[F′μ(U(p))]-1Fμ(X(p))

(26)

X(p+1)=V(p)+{[F′μ(X(p))]-1-

2[F′μ(U(p))]-1}Fμ(V(p))

(27)

X(p)=X(p+1),p=p+1,轉式(20)。

(3)輸出計算結果。

z*=|x(p)|-x(p),w*=|x(p)|+x(p),y*=y(p)

(28)

文獻[16]證明了該方法在求解非線性方程組時具有五階收斂速度；文獻[17]應用該方法求解了線性規劃,首次把該方法應用于求解凸二次規劃。在實際計算過程中,若雅克比矩陣接近奇異,可以采用阻尼牛頓法處理。

4 數值實驗

采用牛頓法求解凸二次規劃的問題,程序采用MATLAB R2009a編寫。設置μ=1×10-3,x(0)=0∈Rn,y(0)=0∈Rm,終止條件ε=1×10-6。

算例1考慮如下凸二次規劃：

(29)

算法經過2次迭代,得z*=(1.129 0, 0.774 2, 0.096 8, 0)T,w*=(0, 0, 0, 1.032 3)T,y*=(0, -1.032 3)T,相應的目標值為-7.161 3。

算例2考慮如下凸二次規劃:

(30)

算法經過2次迭代,得z*=(1, 1, 0)T,w*=(0, 0, 2)T,y*=-2,相應的目標值為-6。

算例3考慮如下凸二次規劃：

(31)

算法經過3次迭代,得z*=(1, 0, 2, 5)T,w*=(0, 2, 0, 0)T,y*=(0, 0)T,相應的目標值為-1。

算例4考慮如下凸二次規劃：

(32)

算法經過3次迭代,得z*=(0, 0, 5)T,w*=(0, 4, 0)T,y*=0,相應的目標值為0。

算例5考慮如下凸二次規劃：

(33)

算法經過2次迭代,得z*=(0, 5, 5)T,w*=(16, 0, 0)T,y*=(-15, 47)T,相應的目標值為195。

算例6考慮如下凸二次規劃：

(34)

算法經過2次迭代,得z*=(6.333 3, 5, 11.666 7, 0)T,w*=(0, 0, 0, 4.666 7)T,y*=(-0.666 7, 4.333 3)T,相應的目標值為-168.166 7。

5 結論

通過把凸二次規劃轉化為非線性方程組,采用光滑逼近函數進行處理,得到光滑非線性方程組,進而利用高階牛頓法進行求解,選用的牛頓法具有五階收斂性；進一步可以構造逼近程度更好的一致光滑函數,采用更高階的牛頓法[18-24]進行求解。