求解絕對(duì)值方程的光滑牛頓算法

張焱嬌

(天津大學(xué)，天津300072)

絕對(duì)值方程一般數(shù)學(xué)表達(dá)式如下:

其中:|x|表示的是對(duì)x的每個(gè)分量取絕對(duì)值.本文中研究的矩陣A和B矩陣為方陣:

絕對(duì)值方程是Rohn在中首次出提出的.目前求解絕對(duì)值方程的算法針對(duì)B=－I，

[1－5]中討論了絕對(duì)值方程(1)(2)和(3)的一些性質(zhì).下面簡(jiǎn)單列舉一下這些性質(zhì):如果絕對(duì)值等式(1)、(2)、(3)是可解的，則絕對(duì)值方程或者有惟一解，或者有多重解[3]，一般的絕對(duì)值規(guī)劃(1)可以轉(zhuǎn)換成一個(gè)線性互補(bǔ)問(wèn)題[4].

Mangasarian提出算法是:將絕對(duì)值等式(3)轉(zhuǎn)化成分段線性的凹極小化問(wèn)題，利用連續(xù)化線性算法來(lái)求解.之后在參考文獻(xiàn)[6]又中提出了一種廣義牛頓法來(lái)求解絕對(duì)值方程(3)，要求A的奇異值嚴(yán)格大于1，此方法的有線性收斂性質(zhì)已經(jīng)得到了證明.之后Oleg Prokopyev在參考文獻(xiàn)[7]中將絕對(duì)值方程(3)問(wèn)題轉(zhuǎn)化為0－1混合線性規(guī)劃來(lái)進(jìn)行求解.本文所做的工作是利用光滑牛頓算法來(lái)求解問(wèn)題(2)和(3)，光滑化的思想[8－9].

本文中證明了在矩陣A的主對(duì)角線上的元素的絕對(duì)值大于矩陣B的主對(duì)角元素的絕對(duì)值，矩陣A、B對(duì)于可交換條件下，算法是適定的，具有全局收斂性質(zhì).根據(jù)參考文獻(xiàn)[10－11]，絕對(duì)值方程的可解性問(wèn)題是一個(gè)NP難的問(wèn)題，所以本文不討論絕對(duì)值方程解的存在性.

現(xiàn)將本文所使用的符號(hào)如下:對(duì)于向量x∈Rn，若x＞0表示的是向量x的每個(gè)分量都大于0;‖x‖表示向量x的歐式范數(shù)，對(duì)于向量y，s∈ Rn對(duì)于函數(shù)F:Rn→Rm，F(xiàn)'(x)表示函數(shù)F在點(diǎn)的Jacobi矩陣.

1 算法設(shè)計(jì)

1.1 NCP 函數(shù)

定義 設(shè)函數(shù)φ(a，b):R×R→R，如果滿足φ(a，b)=0?a≥0，b≥0，且 ab=0 性質(zhì).

1.2 光滑牛頓法

定理1 對(duì)于絕對(duì)值等式Ax+B|x|=b，其中x∈Rn，A∈Rn×n，B∈Rn×n若矩陣(－ A+B)和矩陣B可逆，矩陣A和矩陣B對(duì)于乘法可交換時(shí)，則該問(wèn)題可等價(jià)轉(zhuǎn)換為如下:

(－A+B)s－(－A－B)y=2Bb y≥0，s≥0，yTs=0

證明:見(jiàn)參考文獻(xiàn)[5]

推論1 當(dāng)B=－I時(shí)，求Ax－|x|=b解的解，若矩(A+I)陣可逆，那么(2)等價(jià)轉(zhuǎn)化為求解如下的問(wèn)題:

現(xiàn)在在利用光滑后的F－B函數(shù)，來(lái)構(gòu)造向量值函數(shù):

Ψ(μ，y，s)=(φ(μ，y1，s1)，φ(μ，y2，s2)，…，φ(μ，yn，sn))T

定義的光滑化函數(shù)H:R2n+1→R(2n+1)，如下:其中 z=(μ，y，s)，顯然‖H(z)‖ =0當(dāng)且僅當(dāng) μ=0，且(y，s)是線性互補(bǔ)問(wèn)題(4)的解.只要求出z*使得‖H(z*)‖=0，就可以求出絕對(duì)值方程(2)的解.對(duì)于任意(μ，y，s)∈R++× Rn× Rn，函數(shù) H(z)一定是連續(xù)可微的.且函數(shù)H(z)Jacobi矩陣為:

下面對(duì)一些符號(hào)予以說(shuō)明:

算法1(求解形如(2)的絕對(duì)值等式的光滑牛頓法)

步驟0選擇適當(dāng)?shù)?δ，σ∈(0，1)，μ0＞0 和(y0，s0)∈R2n是初始向量.令 z0:=(μ0，y0，s0)，選擇適當(dāng)?shù)摩拢?使得‖H(z0)‖≤βμ0，令 e0=(1，0，…，0)T，令 k:=0;

步驟1如果‖H(zk)‖=0，停止，zk則是我們所求的最優(yōu)解;否則，跳至步驟2;

步驟2 計(jì)算 Δzk:=(Δμ0，Δy0，Δs0)∈R ×Rn×Rn，其中Δzk滿足以下的牛頓方程:

根據(jù)(5)的具體形式，可以得到:

步驟3令λk是1，δ，δ2中使下式成立的最大值.

步驟4 令 zk+1:zk+λkΔzk取 k:=k+1，調(diào)到步驟1.

步驟5求解方程y－s=2Bx

1.3 算法的性質(zhì)

在參考文獻(xiàn)[8]中，有如下結(jié)論成立:

命題1 設(shè){zk=(μk，yk，sk)是由算法 1 產(chǎn)生的迭代序列，則有以下的結(jié)論成立:

迭代過(guò)程中，對(duì)任意迭代點(diǎn)都有μk＞0;序列{‖H(zk)‖}是單調(diào)下降序列;算法迭代的過(guò)程中序列{μk}是單調(diào)下降的;設(shè)β＞1是算法中給定的常數(shù)，定義集合:N(β):={z∈R+×Rn×Rn:‖H(z)‖≤βμ}則我們?cè)诘^(guò)程中，對(duì)于任意的k∈{1，2，…，n}，zk∈N(β);

定理2若果矩陣(A+B)，(－A+B)和矩陣B可逆，矩陣A和B矩陣對(duì)于乘法可交換，并且矩陣A的主對(duì)角線上的元素的絕對(duì)值嚴(yán)格大于矩陣B的主對(duì)角線上元素的絕對(duì)值，則算法1中的步驟2的牛頓方程一定是可解的，說(shuō)明此算法一定是適定的.

證明:先介紹一些證明中涉及的符號(hào):

cii:=min{－ aii－ bii，aii+bii}，

c=(c11，c22，cnn)C:=diag(c)

dii:=max{－ aii－ bii，aii+bii}，

d=(d11，d22，knn)D:=diag(d)，i∈{1，2，…，n}

由假設(shè)矩陣A的對(duì)角線上的元素的絕對(duì)值大于對(duì)應(yīng)B的位置元素的絕對(duì)值，并且下式成立:

|aii－1|bii＜1，i∈{1，2，…，n}

由參考文獻(xiàn)[13]，推論 3.2區(qū)間[C，D]是正則的.即任意的，m，n∈[0，1]，mC+nD 是可逆的.

Ψs(A+B)－Ψy(－A+B)是可逆的

?Ψs(A+B)－Ψy(A+B)－1(A+B)(－A+B)是可逆的

?Ψs(A+B)－Ψy(A+B)－1(A+B)(－A+B)是可逆的，因?yàn)锳B=BA

?Ψs－Ψy(A+B)－1(－A+B)是可逆的，因?yàn)榫仃囀?A+B)可逆的.

那么算法1的步驟2的牛頓方程也一定有解的.

推論2當(dāng)矩陣B=－I時(shí)，絕對(duì)值方程(3)有以下的結(jié)論成立:如果矩陣(I－A)可逆，并且矩陣的主對(duì)角線上的元素同時(shí)大于1或者小于－1，則算法中的牛頓方程是可解的.

2 算法的性質(zhì)

本節(jié)證明在矩陣AB滿足適當(dāng)關(guān)系的條件下，算法1有全局收斂性.

定理3矩陣(I－A)，(I+A)是可逆的，并且矩陣A的主對(duì)角線上的元素同時(shí)嚴(yán)格大于1或者嚴(yán)格小于－1時(shí)，則通過(guò)算法1產(chǎn)生的{zk:=(μk，yk，sk)}是有界序列.

證明 由命題1的結(jié)論可知序列{μk}是單調(diào)下降有下界的序列{zk}:={μk，wk}，為了說(shuō)明序列的有界性{wk}，只需在說(shuō)明序列是有界序列即可，其中{wk=(yk，sk)}.利用反證法證明序列{wk}是是有界序列{wk}，假設(shè)序列無(wú)界，則可知其部分序列{yk}，{sk}至少有一個(gè)是無(wú)界序列.則無(wú)界序列，必須滿足以下三種情況之一:1)序列{yk}無(wú)界，序列{sk}有界;2)序列{yk}有界，序列{sk}無(wú)界;3)序列{yk}，{sk}無(wú)界;

1)序列{yk}無(wú)界，序列{sk}有界的

因?yàn)樾蛄衶yk}無(wú)界，則存在指標(biāo) i0∈{1，2，…，n}，使得|(yk)|→∞，當(dāng) k→∞時(shí)，

|((I+A)sk+(I－A)yk+2b)i0|→∞，k→∞，

此時(shí)就一定會(huì)出現(xiàn)‖H(zk)‖→∞，當(dāng)k→∞，這與命題2.1中序列‖H(zk)‖是單調(diào)下降有界的序列產(chǎn)生矛盾.所以情況(1)是不可能出現(xiàn)的.

2)序列{yk}有界，序列{sk}無(wú)界的(同(1))

3)序列{yk}，{sk}同時(shí)是無(wú)界的

若實(shí)數(shù)對(duì){i0，j0}其中 i0∈{1，2，…，n}，j0∈{1，2，…，n}且滿足i0≠j0則此時(shí)可以歸結(jié)為上述的(a)、(b)兩種情況之一，會(huì)出現(xiàn)矛盾，無(wú)界序列無(wú)法產(chǎn)生.下面討論另外一種情況:即為對(duì)于任意t0∈{1，2，…n}，只要|(yk)t0|→∞ 則必有|(sk)t0|→∞并且其逆命題也成立.下面在進(jìn)一步討論，當(dāng)k→∞時(shí):

若(yk)t0→ +∞，(sk)t0→ +∞ 時(shí)，((I+A)sk+(I－A)yk+2b)t0→+∞，

若(yk)t0→ +∞，(sk)t0→ +∞ 時(shí)，((I+A)sk+(I－A)yk+2b)t0→+∞

若(yk)t0→+∞，(sk)t0→－∞時(shí)，有

顯然當(dāng) k→∞，有|φ(μ，(yk)t0，(sk)t0|→∞，則此時(shí)同樣可以的得到函數(shù)H(zk)的范數(shù)‖H(zk)‖→∞，與算法性質(zhì)矛盾.綜上所述，假序列{zk:=μk，yk，sk)}也有界的.

定理4 假定矩陣(I－A)，(I+A)是可逆的，矩陣A的對(duì)角線上的元素同時(shí)大于1或者同時(shí)小于－1，則產(chǎn)生的序列{zk}有以下性質(zhì)

證明方法類似于參考文獻(xiàn)[13]

3 數(shù)值計(jì)算結(jié)果與結(jié)論

數(shù)值實(shí)驗(yàn)所用電腦規(guī)格為:CPU是2.93 GHz，內(nèi)存2.00 GHz.數(shù)值實(shí)驗(yàn)的具體絕對(duì)值方程產(chǎn)生過(guò)程如下:產(chǎn)生一個(gè)矩陣保證其每一個(gè)對(duì)角元素都同時(shí)大于1或者同時(shí)小于－1;隨機(jī)選取向量x∈Rn，其中的每個(gè)分量符合[－10，10]的均勻分布;計(jì)算出b=Ax－|x|.

本文將對(duì)n等于500、800、1000的情況分別進(jìn)行數(shù)值實(shí)驗(yàn)，每組實(shí)驗(yàn)中，連續(xù)隨機(jī)產(chǎn)生50個(gè)絕對(duì)值方程.在數(shù)值實(shí)驗(yàn)的程序中‖Ax－|x|－b‖＜10－6作為算法的終止條件;循環(huán)的最高次數(shù)定為40次;當(dāng)算法終止是若‖Ax－ |x|－b‖ ＜10－6，則認(rèn)為失敗.

相關(guān)參數(shù)值為 σ =0.8，σ =0.0001，μ0=0.0001;x0為隨機(jī)產(chǎn)生的維向量x0i，其分量符合[－10，10]之間的均勻分布.

表1 數(shù)值結(jié)果(矩陣的主對(duì)角線上的元素嚴(yán)格大于1或者嚴(yán)格小于－1)

4 結(jié)語(yǔ)

本文提出的光滑化牛頓算法，分別對(duì)形如Ax+B|x|=b和Ax－|x|=b的絕對(duì)值方程組進(jìn)行求解，并且證明了在一定的條件下此算法是全局收斂的.在此之前，要求的是在區(qū)間[I－A，I+A]是正則的，對(duì)矩陣A的奇異值等性質(zhì)要求較高.

在理論證明中，條件“當(dāng)矩陣(B+A)和矩陣(B－A)可逆，矩陣AB對(duì)乘法滿足交換律，以及矩陣A的主對(duì)角線上的元素的絕對(duì)值大于矩陣B主對(duì)角線上的絕對(duì)值時(shí)，證明了算法是適定的，并且當(dāng)矩B=－I陣時(shí)，更證明了該算法的全局收斂性，”在未來(lái)的研究中，可以嘗試證明矩陣B的一般情形是算法的全局收斂性.

參考文獻(xiàn):

[1]MANGASARIAN O L.Knapsack feasibility as an absolute value equation solvable by successive linear programming[J].Optimization Letters，2009，3(2):161 －170.

[2]OLEG P.On equivalent reformulations for absolute value equations[J].Computational Optimization and Applications，2009，44(3):363－372.

[3]FUKUSHIMA M.The primal Douglas－Rachford splitting algorithm for a class of monotone mappings with application to the traffic equilibrium problem [J].Mathematical Programming，1996，72(1):1－15.

[4]NAGURNEY A，DONG J，ZHANG D.A supply chain network equilibrium model[J].Transportation Research，2002，Part E 38(5):281－303

[5]PANG TRINKLE J C.Complementarity formulations and existence of solutions of multi－rigid－body contact problems with Coulomb friction[J].Mathematical Programming，1996，73(2):199－226.

[6]CHEN CH，MANGASARIAN O L.Smoothing methods for convex inequalities and linear complementarity problems[J].Mathematical Programming，1995，71:51－69.

[7]MANGASARIAN O L，MEYER R.Absolute Value Equations[J].Linear Algebra and Its Applications，2006，419:359 －367.

[8]CHEN C， MANGASARIAN O L.A class of smoothing functions for nonlinear and mixed complementarity problems[J].Computational Optimization and Applications，1996，5:97－138.

[9]QI L Q，SUN D，Smoothing functions and smoothing Newton method for complementarity and variational inequality problems[J].Jouranal of Optimization Theory and Applications，2002，113 －121.

[10]ROHN J.A theorem of the alternatives for the equation Ax+B|x|=b[J].Linear Multilinear Algebra，2004，52(6):421－426.

[11]MANGASARIAN O L .Absolute value programming[J].Computational Optimization and Applications，2007，36:43 －53.

[12]MANGASARIAN O L.Absolute value equation solution via concave minimization[J].Optimization Letters，2007，1(1):3－8.

[13]高竹峰.求解絕對(duì)值方程組的牛頓算法[D].天津:天津大學(xué)，2009.