求解非負(fù)約束優(yōu)化問題的可行LS共軛梯度法*

劉陶文,周 蓉

( 湖南大學(xué) 數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院,湖南 長沙 410082 )

非負(fù)約束優(yōu)化問題廣泛存在于許多學(xué)科及工程應(yīng)用領(lǐng)域中,如:非線性回歸分析、金融投資、大地測(cè)量、衛(wèi)星導(dǎo)航等,并且有關(guān)的數(shù)據(jù)處理是非常龐大的,呈現(xiàn)的問題通常是大規(guī)模的.因此,研究求解這些問題的高效算法具有重要的現(xiàn)實(shí)意義.非負(fù)約束優(yōu)化問題的一般形式為

(1)

這里函數(shù)f:Rn→R是連續(xù)可微的,D={x∈Rn:x≥0}是可行域.

問題(1)的一階最優(yōu)性條件 ( KKT條件 )可以表示為

(2)

這里g(x)=▽f(x)表示函數(shù)f在點(diǎn)x處的梯度,gi(x)表示g(x)的第i個(gè)分量．滿足條件(2)的點(diǎn)x被稱為問題(1)的一個(gè)KKT點(diǎn)．

與問題(1)相關(guān)的是下面的有界約束問題：

(3)

這里l和u是Rn中的有界向量．許多研究集中于求解大規(guī)模有界約束問題(3),人們提出了如譜梯度投影法、有限記憶擬牛頓法和共軛梯度法等有效算法[1,4,6]．盡管問題(1)可以看作問題(3)的特殊情形,但不知這些方法能否經(jīng)過適當(dāng)?shù)男薷谋挥脕砬蠼夥秦?fù)約束問題．

眾所周知，共軛梯度法及其修正形式是求解大規(guī)模無約束問題minf(x),x∈Rn的有效方法[2,7-8]．因此,人們?cè)O(shè)想用共軛梯度法的思想求解約束問題.文獻(xiàn)[6]研究了用共軛梯度的思想求解有界約束問題(3),通過求解一個(gè)約束子問題來計(jì)算搜索方向,作者證明了在 Wolfe型線性搜索下所提出的方法是收斂的．然而已有研究表明在該研究上存在許多的困難[5].

在本文中，我們研究用共軛梯度型方法求解非負(fù)約束問題 (1).結(jié)合已有求解有界約束問題的子空間思想[4]和求解無約束問題的Liu-Storey (LS)共軛梯度法[3]，提出了一個(gè)求解問題(1)的可行LS共軛梯度法,與文獻(xiàn)[6]的方法比較,本文提出的方法的優(yōu)點(diǎn)是無需求解子問題,并且在Armijo搜索下具有全局收斂性,節(jié)省大量計(jì)算．

1 計(jì)算下降可行方向

給定x∈D及常數(shù)ε>0,首先定義如下的索引集合：

B(x)={i:xi>ε}，

A1(x)={i:xi=0,gi(x)≥0}∪{i:0

A2(x)={i:0≤xi≤ε,gi(x)<0}，

A3(x)={i:0

(4)

我們稱集合A(x)=A1(x)∪A2(x)∪A3(x)為x處的有效集,而B(x)稱為非有效集．

設(shè)xk∈D是問題(1)的一個(gè)近似解,為了計(jì)算在該點(diǎn)處的下降可行搜索方向dk∈Rn,采用如下策略：由于A1(xk)中的變量滿足條件(2),這些變量的取值保持不變; 對(duì)于A3(xk)中的變量,由于-gi(xk)指向可行域的邊界,對(duì)-gi(xk)進(jìn)行截?cái)嗵幚慝@得搜索方向；對(duì)于B(xk)∪A2(xk)中的變量先計(jì)算某種共軛梯度型方向,然后進(jìn)行適當(dāng)?shù)慕財(cái)嗵幚懋a(chǎn)生搜索方向.dk的具體計(jì)算方法如下：

設(shè)矩陣Pk由列{ei:i∈B(xk)∪A2(xk)}構(gòu)成,矩陣Qk由列{ek:i∈A3(xk)}構(gòu)成,其中ei表示n階單位矩陣的第i列.首先計(jì)算共軛梯度方向如下:

(5)

這里gk=g(xk),yk=gk-gk-1,并且

(6)

容易驗(yàn)證

令

(7)

(8)

(9)

因此,由dk的定義(7)可得

(10)

證畢

證由KKT條件(2)以及(10)即可驗(yàn)證定理結(jié)論成立.

證畢

2 算法結(jié)構(gòu)及其收斂性分析

現(xiàn)在提出求解問題(1)的可行LS共軛梯度算法如下:

算法1(可行LS共軛梯度法)

步 0:取初始點(diǎn)x0∈D,常數(shù)ρ,σ∈(0,1),ε>0.令k=0.

加拿大林間貓頭鷹的鳴聲“呵呵——呵——呵——”，極像老翁冷冷的笑聲，與我記憶中江南的貓頭鷹的叫法“骨碌碌碌……”圓滾而陰森的鳴音，大異其趣。鴟鸮的鳴法會(huì)超過人類的方言上萬種嗎？鳥啼人語，都是天趣。

f(xk+αdk)≤f(xk)-σ‖αdk‖2

(11)

的最大值αk.令xk+1=xk+αkdk.

步4:令k:=k+1,然后轉(zhuǎn)步1.

在算法1的步3中,使用了Armijo 型線性搜索(11)計(jì)算步長αk,這比文獻(xiàn)[6]中使用的Wolfe 型搜索更容易實(shí)現(xiàn),而且節(jié)省大量投影計(jì)算.顯然由算法1產(chǎn)生的序列{xk}?D．

下面分析算法1的全局收斂性.首先作如下的假設(shè)：

假設(shè)A由算法1產(chǎn)生的序列{xk}在D內(nèi)是有界的,函數(shù)f(x)在D內(nèi)有下界并且其梯度g(x)在可行域D內(nèi)是Lipschitz連續(xù)的,即存在常數(shù)L>0滿足

‖g(x)-g(y)‖ ≤L‖x-y‖,?x,y∈D.

注當(dāng)dk≠0,dk是f(x)在xk處的一個(gè)下降可行方向,從而由假設(shè)A可知序列{f(xk)}單調(diào)下降且有下界,從而得到

(12)

(13)

證根據(jù)索引分類(4)分4種情形進(jìn)行討論．

情形1：當(dāng)i∈B(xk)時(shí),(xk)i>ε,故

情形2：當(dāng)i∈A1(xk)時(shí),(xk)i+γ(dk)i=0,?γ≥0.

情形3：當(dāng)i∈A2(xk)時(shí),

情形4：當(dāng)i∈A3(xk)時(shí),gi(xk)≥0,即有

以上的討論表明(13)成立．

證畢

引理3如果存在某個(gè)常數(shù)δ>0使得

那么序列{dk}是有界的,即存在常數(shù)M>0使得‖dk‖∞≤M.

證在假設(shè)A的條件下,{xk},{gk}是有界的,即存在常數(shù)M1>0,使得

‖xk‖∞≤M1，‖gk‖∞≤M1,?k≥0.

(14)

另一方面,由(12),必存在常數(shù)μ∈(0,1)及正整數(shù)k1,使得

2δ-1M1L‖αk-1dk-1‖ ≤μ,?k≥k1.

因此,?k≥k1,進(jìn)一步有

所以,對(duì)任意的k≥0,有

證畢

利用引理3,可以得到下面的全局收斂性定理.

定理2設(shè)序列{xk}由算法1產(chǎn)生.如果假設(shè)A 成立,則有

(15)

下面對(duì)步長αk分兩種情形討論．

另一方面,由泰勒中值定理及g(x)的 Lipschitz 連續(xù)性容易驗(yàn)證

由上面的兩個(gè)不等式可推出

證畢

由定理1及定理2知,必存在迭代序列{xk}的某一個(gè)極限點(diǎn)x*是問題 (1)的KKT點(diǎn)．

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

在這一節(jié)，將所提出的算法應(yīng)用于大地測(cè)量中數(shù)據(jù)處理問題[9].對(duì)一理想邊坡因地質(zhì)斷層構(gòu)成一可能的滑體,按地質(zhì)學(xué)知識(shí),滑體只能沿底盤向東北方向偏下移動(dòng).選定三個(gè)基點(diǎn),其坐標(biāo)分別為A(500.00,0,100.00),B(0.00,-500.00,150.00),C(500.00,500.00,200.00).在滑體上,取監(jiān)測(cè)點(diǎn)P(0,0,100.00),分別在基點(diǎn)A，B，C處用邊前方交會(huì)監(jiān)測(cè)P點(diǎn)的位移(x,y,-z),AP，BP和CP為3個(gè)觀測(cè)邊,且有先驗(yàn)信息x,y,z>0.設(shè)觀測(cè)向量為L∈R3,則有相應(yīng)的誤差方程:

L+V=(|AP|,|BP|,|CP|)T,

其中V為觀察噪聲向量.設(shè)A，B，C3個(gè)邊觀測(cè)的一次觀測(cè)值分別為500.04,502.52,714.13.試估計(jì)P點(diǎn)的位移(x,y,-z).

當(dāng)P點(diǎn)有位移(x,y,-z)時(shí),觀測(cè)向量L的誤差方程為:

文獻(xiàn)[9]將誤差方程線性化后,利用求解線性互補(bǔ)問題的Lemke算法求解得P點(diǎn)位移的最小二乘估計(jì).由于線性化產(chǎn)生的誤差及Lemke算法的計(jì)算復(fù)雜性,這種估計(jì)的精度是有限的.

記函數(shù).

則P點(diǎn)位移的估計(jì)問題等價(jià)于下面的帶非負(fù)約束的非線性最小二乘問題:

(16)

[1]BIRGIN E G,MARTINEZ J M,RAYDAN M.Non-monotone spectral projection gradient methods on convex sets [J].SIAM Journal on Optimization,2000,10(1):1196-1211.

[2]DAI Y,YUAN Y.Nonlinear conjugate gradient methods [M].Shanghai:Shanghai Science and Technology Publisher,2000.

[3]LIU Y,STOREY C.Efficient generated conjugate gradient algorithms Part 1:Theory [J].Journal of Optimization Theory and Applications,1991,69(1):129-137.

[4]NI Q,YUAN Y.A subspace limited memory quasi-Newton algorithm for large-scale nonlinear bound constrained optimization [J].Mathematics of Computation,1997,66(220):1509-1520.

[5]NOCEDAL J.Conjugate gradient methods and nonlinear optimization,In:Linear and nonlinear conjugate gradient-Related method [M].Philadelphia,SIAM Press,1996,9-23.

[6]PYTLAK R.An efficient algorithmfor large-scale nonlinear programming problems with simple bounds on the variables [J].SIAM Journal on Optimization,1998,8(2):532-560.

[7]ZHANG L,ZHOU W J,LI D H.A descent modified Polak-Rebière-Polyak conjugate gradient method and its global convergence [J].IMA Journal of Numerical Analysis,2006,26(4):629-640.

[8]ZHANG L,JIAN S Y.Further studies on the Wei-Yao-Liu nonlinear conjugate gradient method [J].Applied Mathematics and Computation,2013,219 (14):7616-7621.

[9]宋迎春,朱建軍,羅德仁，等.附非負(fù)約束平差模型的最小二乘估計(jì) [J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2008,33(9):907-909.

SONG Ying-cun,ZHU Jian-jun,LUO De-ren,etal.Least-squares estimation of nonnegative constrained adjustment model [J].Geomatics and Information Science of Wuhan University:Information and Science,2008,33(9):907-909 ．(In Chinese)