李丹丹, 李遠(yuǎn)飛, 王松華

(1. 廣州華商學(xué)院 應(yīng)用數(shù)學(xué)系, 廣州 511300; 2. 百色學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 廣西 百色 533000)

0 引 言

非線性方程組廣泛應(yīng)用于通信、 化學(xué)、 電力系統(tǒng)、 壓縮感知等領(lǐng)域[1-3], 因此建立高效求解非線性方程組的數(shù)值方法具有重要意義. 本文主要考慮如下非線性方程組：

H(z)=0,z∈n,

(1)

其中H:n→n為連續(xù)單調(diào)函數(shù), 即對(duì)于任意的z1,z2∈n, 有

(H(z1)-H(z2))T(z1-z2)≥0.

(2)

目前, 求解非線性方程組的方法主要有牛頓法、 擬牛頓法、 共軛梯度法等[4-5], 其中牛頓法和擬牛頓法由于局部收斂較快等優(yōu)勢(shì), 廣泛應(yīng)用于求解中小規(guī)模的非線性方程組中, 但在運(yùn)算過(guò)程中每次迭代都需要計(jì)算一個(gè)與Jacobi矩陣相關(guān)的線性方程組, 運(yùn)算復(fù)雜性較高, 限制了其在求解大規(guī)模非線性方程組問(wèn)題中的應(yīng)用.共軛梯度法由于在求解過(guò)程中無(wú)需計(jì)算方程組的Jacobi矩陣、 編程簡(jiǎn)單且存儲(chǔ)空間要求較小等特點(diǎn)而成為求解大規(guī)模優(yōu)化問(wèn)題的主要方法之一, 目前已成功地應(yīng)用于求解大規(guī)模非線性方程問(wèn)題和壓縮感知問(wèn)題中[6-8].

不失一般性, 共軛梯度法的迭代式為zk+1=zk+tkdk, 其中:zk為當(dāng)前迭代點(diǎn);zk+1為下一個(gè)迭代點(diǎn);tk為某種線搜索產(chǎn)生的步長(zhǎng);dk為某種方法產(chǎn)生的搜索方向, 其迭代形式為

式中βk為共軛參數(shù),H(zk)簡(jiǎn)記為Hk.顯然, 共軛梯度法的優(yōu)劣受給定線搜索方法和搜索方向的影響, 而搜索方向的好壞取決于共軛參數(shù), 其決定了算法的收斂速度.

針對(duì)單調(diào)非線性優(yōu)化問(wèn)題, 文獻(xiàn)[9]采用一種高效的線搜索技術(shù), 令步長(zhǎng)tk=srmk,mk為滿足下列不等式的最小非負(fù)整數(shù)：

-H(zk+tkdk)Tdk≥βtk‖dk‖2,

(3)

其中s>0, 00.

經(jīng)典Hestenes-Stiefel(HS)方法因自動(dòng)滿足共軛條件、 數(shù)值結(jié)果較好等優(yōu)勢(shì)而備受關(guān)注, 但研究表明, 該方法收斂性質(zhì)一般. 為改善HS算法的收斂性質(zhì), 文獻(xiàn)[10]在經(jīng)典HS方法的基礎(chǔ)上, 通過(guò)添加一個(gè)擾動(dòng)項(xiàng), 給出了一種三項(xiàng)搜索方向, 其表達(dá)式為

(4)

Rk={z∈n|H(uk)T(z-uk)=0},

將最優(yōu)解z*與當(dāng)前迭代點(diǎn)zk嚴(yán)格分離, 下一個(gè)新的迭代點(diǎn)zk+1可由下式計(jì)算得到：

(5)

基于上述分析, 本文在式(4)的基礎(chǔ)上, 通過(guò)構(gòu)建一個(gè)具有良好性質(zhì)的新型搜索方向, 利用式(3)的線搜索技術(shù)和式(5)的超平面投影方法, 提出一種無(wú)導(dǎo)數(shù)修正三項(xiàng)HS共軛梯度投影方法, 分析新算法的全局收斂性, 給出數(shù)值實(shí)驗(yàn)結(jié)果, 并將其應(yīng)用于信號(hào)恢復(fù)問(wèn)題中.

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

(6)

因此, 式(4)可改寫(xiě)為

(7)

(8)

其中δ>0為給定的常數(shù).式(8)所構(gòu)造的搜索方向不僅保留了文獻(xiàn)[10]的充分下降性, 而且克服了文獻(xiàn)[10]中定義的不適定性等問(wèn)題.

通過(guò)上述討論, 可建立一種修正三項(xiàng)HS共軛梯度投影算法, 用于求解非線性單調(diào)方程組(1), 算法描述如下.

算法1MHS(Modified Hestenes-Stiefel)算法.

步驟1) 給定初始點(diǎn)z0∈n(k=0), 運(yùn)行參數(shù)s>0, 00,δ>0；

步驟2) 計(jì)算函數(shù)值H(zk), 如果‖H(zk)‖≤ε, 則算法停止, 否則轉(zhuǎn)步驟3)；

步驟3) 利用式(6)和式(8)計(jì)算搜索方向dk；

步驟4) 利用式(3)的線搜索方法計(jì)算步長(zhǎng)tk, 令uk=zk+tkdk；

步驟5) 計(jì)算函數(shù)值H(uk), 如果‖H(uk)‖≤ε, 則算法停止, 否則由式(5)計(jì)算出新的迭代點(diǎn)zk+1, 令k∶=k+1, 轉(zhuǎn)步驟2).

2 全局收斂性分析

下面分析MHS算法的全局收斂性質(zhì).

假設(shè):

(H1) 非線性單調(diào)方程組(1)的解集非空；

(H2) 函數(shù)H(z)在n上是Lipschitz連續(xù)的, 即存在常數(shù)a>0, 使得

‖H(z1)-H(z2)‖≤a‖z1-z2‖, ?z1,z2∈n.

(9)

引理1若搜索方向dk由式(6)和式(8)產(chǎn)生, 則搜索方向dk滿足充分下降性, 即對(duì)于任意的k, 均有

(10)

于是式(10)得證.

利用Cauchy-Schwatz不等式, 由式(10)易推出

‖Hk‖≤‖dk‖.

(11)

下面證明MHS算法的適定性并給出步長(zhǎng)tk的下界.

引理2如果假設(shè)(H1),(H2)成立, 則存在一個(gè)步長(zhǎng)tk使得式(3)成立.

證明： 類似于文獻(xiàn)[13]中引理3.1可證, 故略.

引理3如果假設(shè)(H1),(H2)成立, 且MHS算法產(chǎn)生序列{zk}和{uk}, 則有

聯(lián)合式(10)和式(9), 得

證明： 先證序列{zk}和{uk}有界.假設(shè)方程組(1)的最優(yōu)解為z*, 即H(z*)=0, 由式(2)得(H(uk)-H(z*))T(uk-z*)≥0, 可化簡(jiǎn)為H(uk)T(zk-z*+uk-zk)≥0, 于是有

H(uk)T(zk-z*)≥H(uk)T(zk-uk).

(12)

此外, 由式(3)的線搜索方法及uk的定義得

(13)

觀察式(14)可知, {‖zk-z*‖}是關(guān)于k的單調(diào)不增序列且收斂, 因此序列{zk}有界.由式(14)還可得：

‖zk+1-z*‖2≤‖zk-z*‖2,

(15)

利用遞推關(guān)系可得

‖zk-z*‖2≤‖zk-1-z*‖2≤…≤‖z0-z*‖2,

結(jié)合假設(shè)(H1),(H2)可得

‖H(zk)‖=‖H(zk)-H(z*)‖≤a‖z0-z*‖.

令c=a‖z0-z*‖, 則序列{H(zk)}有界, 即對(duì)于任意的k, 均有

‖H(zk)‖≤c.

(16)

由式(13)和Cauchy-Schwatz不等式得

再聯(lián)合式(16)和序列{zk}的有界性, 可知序列{uk}也有界.

因?yàn)樾蛄衶uk}有界, 所以{‖uk-z*‖}也有界, 故存在一個(gè)正常數(shù)e, 使得‖uk-z*‖≤e成立, 再結(jié)合假設(shè)(H1),(H2)得

‖H(uk)‖=‖H(uk)-H(z*)‖≤a‖uk-z*‖≤ae,

聯(lián)合式(14)得

從而有

下面證明MHS算法的全局收斂性質(zhì).

證明： 用反證法.假設(shè)結(jié)論不成立, 則存在κ>0, 使得對(duì)任意的k, 均有

‖H(zk)‖≥κ.

(17)

結(jié)合式(11)得

κ≤‖Hk‖≤‖dk‖.

(18)

另一方面, 由引理4可知{zk}有界, 故存在正常數(shù)e1, 使得‖zk‖≤e1, 再結(jié)合wk-1的定義及假設(shè)(H1),(H2)得

‖wk-1‖=‖Hk-Hk-1‖≤a‖zk-zk-1‖≤a(‖zk‖+‖zk-1‖)≤2ae1.

(19)

(20)

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

3.1 非線性單調(diào)方程組

為考察MHS算法的性能, 下面將其與MCG(Modified three-term CG method)算法[13]和MPRP(Modified Polak-Ribière-Polyak)算法[11]進(jìn)行數(shù)值效果對(duì)比. 算法參數(shù)設(shè)置為s=1,r=0.5,β=0.01,δ=1.61. 測(cè)試計(jì)算機(jī)環(huán)境為Windows10(64 bite), RAM 8 GB, CPU 3.60 GHz, 用MATLAB(2014a)軟件進(jìn)行編程. 程序終止準(zhǔn)則為‖Hk‖≤10-5, 迭代次數(shù)上限為3 000, 維數(shù)為4 500,12 000,24 000,30 000,45 000.測(cè)試問(wèn)題選自文獻(xiàn)[14], 函數(shù)名稱和初始點(diǎn)列于表1.實(shí)驗(yàn)結(jié)果列于表2, 其中t為程序運(yùn)行時(shí)間, “—”表示迭代次數(shù)已超過(guò)3 000.

表1 測(cè)試問(wèn)題

表2 實(shí)驗(yàn)結(jié)果

續(xù)表2

由表2可見(jiàn)： MHS算法在規(guī)定的迭代次數(shù)內(nèi)均能成功求解這10個(gè)函數(shù), 而MCG算法和MPRP算法不能很好地求解第5個(gè)函數(shù)； MHS算法比MCG算法和MPRP算法在求解同類問(wèn)題時(shí)迭代次數(shù)更少. 為更直觀展示3種不同算法的優(yōu)劣性, 采用文獻(xiàn)[15]的性能曲線評(píng)價(jià)上述3種算法的性能. 圖1～圖3分別為不同算法的迭代次數(shù)、 函數(shù)計(jì)算次數(shù)和CPU運(yùn)行時(shí)間性能比較結(jié)果, 其中曲線越靠上, 表示算法越穩(wěn)定、 有效. 圖1～圖3中橫坐標(biāo)t表示給定的數(shù)值比率， 縱坐標(biāo)計(jì)算公式如下：

其中P表示測(cè)試集, |P|表示測(cè)試集的問(wèn)題個(gè)數(shù),V表示算法集合,rp,s表示某種指標(biāo)(如運(yùn)行時(shí)間、 迭代次數(shù)和函數(shù)計(jì)算次數(shù)等).

圖1 不同算法的迭代次數(shù)性能比較Fig.1 Performance comparison of number of iterations by different algorithms

圖2 不同算法的函數(shù)計(jì)算次數(shù)性能比較Fig.2 Performance comparison of calculation times of functions by different algorithms

圖3 不同算法的運(yùn)行時(shí)間性能比較Fig.3 Performance comparison of running time by different algorithms

由圖1～圖3可見(jiàn), MHS算法總體上比MCG算法和MPRP算法性能更優(yōu), 以更小的迭代次數(shù)和更短的運(yùn)行時(shí)間求解了非線性單調(diào)方程組問(wèn)題.

3.2 信號(hào)恢復(fù)問(wèn)題

下面用MHS算法求解壓縮感知的信號(hào)恢復(fù)問(wèn)題[16], 并與MPRP算法進(jìn)行比較, 以驗(yàn)證本文算法的性能.

表3為兩種算法在求解信號(hào)恢復(fù)問(wèn)題時(shí)的運(yùn)行時(shí)間(t/s)、 迭代次數(shù)及MSE結(jié)果. 圖4為原始信號(hào)(A)、 觀測(cè)信號(hào)(B)、 用MHS算法(C)和MPRP算法(D)恢復(fù)的信號(hào)對(duì)比結(jié)果. 由表3和圖4可見(jiàn), MHS算法和MPRP算法均能有效地從觀察信號(hào)中恢復(fù)出原始信號(hào), 但MHS算法的求解效率更高.

表3 兩種算法對(duì)稀疏信號(hào)的恢復(fù)結(jié)果

(A) 原始信號(hào)(n=4 096, 非零個(gè)數(shù)為128); (B) 觀察信號(hào); (C) MHS算法(MSE=1.317 11×10-5, 迭代次數(shù)為619, 運(yùn)行時(shí)間為68.984 4 s); (D) MPRP算法(MSE=1.190 12×10-5, 迭代次數(shù)為654, 運(yùn)行時(shí)間為78.937 5 s).圖4 兩種算法對(duì)稀疏信號(hào)的恢復(fù)對(duì)比Fig.4 Comparison of sparse signal recovery by two algorithms

實(shí)驗(yàn)結(jié)果表明, 本文提出的新算法不僅具有很好的理論性質(zhì), 且數(shù)值效果較好, 無(wú)論是在求解大規(guī)模非線性方程組問(wèn)題還是在稀疏信號(hào)恢復(fù)問(wèn)題上均有效.