一種特殊的下降算法

摘 要：求解無約束優(yōu)化問題，常用的方法有下降算法，牛頓法，共軛梯度法等。當(dāng)目標函數(shù)為幾個光滑函數(shù)的和時，一些學(xué)者提出并研究了增量梯度算法。其基本思想是循環(huán)選取單個函數(shù)的負梯度作為迭代方向。增量梯度算法的迭代方向不一定是下降方向，所以不能用下降算法的一維搜索確定步長，因為受限于步長的選擇，收斂效率不高。本文結(jié)合了下降算法和增量梯度算法的思想，提出了分裂梯度法。簡單的說，分裂梯度法循環(huán)考慮單個函數(shù)的負梯度方向，如果這一方向是下降方向，則選擇這一方向為迭代方向;否則選取函數(shù)的負梯度方向為迭代方向。最后通過數(shù)值實驗與最速下降算法、隨機下降算法以及增量梯度算法進行對比，結(jié)果表明對于某些優(yōu)化問題，采用分裂梯度法更有效。

關(guān)鍵詞：無約束優(yōu)化;下降算法;增量梯度法;分裂梯度法;Armijo步長規(guī)則

中圖分類號：O224

文獻標識碼： A

本文研究的是優(yōu)化問題

minx∈Rnf（x）∶=∑i∈Ifi（x）。（1）

其中fi∶Rn→R為一階連續(xù)可微函數(shù)， i∈I∶={1，2，…，m}（m為正整數(shù)）。該問題是一個無約束優(yōu)化問題。它在最優(yōu)控制、機器學(xué)習(xí)、數(shù)據(jù)挖掘等許多實際問題中也有著廣泛應(yīng)用，例如求解最小二乘問題以及無線傳感器網(wǎng)絡(luò)中的分布式優(yōu)化問題、神經(jīng)網(wǎng)絡(luò)訓(xùn)練問題等[1-7]。對于此類問題，我們可采用最速下降算法[8]、下降算法[9]、共軛梯度法[8]和擬牛頓法[8]等。但運用此類算法時，每次迭代都需要計算整個函數(shù)的梯度SymbolQC@

f（即需計算每一個fi的梯度SymbolQC@

fi）?？紤]到目標函數(shù)f的特殊性（由多個函數(shù)之和的構(gòu)成形式），一些學(xué)者提出了增量梯度算法[10-11]來求解問題（1）。我們知道下降方向的選擇對算法的迭代次數(shù)有著至關(guān)重要的作用。（i）最速下降算法是選擇迭代方向為dk=-SymbolQC@

f（xk）;（ii）下降算法則是選擇迭代方向dk滿足dkTSymbolQC@

f（xk）lt;0;（iii）增量梯度算法是循環(huán)選擇單個函數(shù)fi的負梯度方向為迭代方向。（i）是1847年由著名數(shù)學(xué)家Cauchy給出的，它是解析法中最古老的一種，其他解析方法或是它的變形，或是受它的啟發(fā)而得到的，因此它是最優(yōu)化方法的基礎(chǔ)，在文[8]定理313中給出了此算法收斂性定理。由于精確一維搜索準則下的最速下降算法在相鄰兩次迭代過程中的迭代方向相互垂直，因而整個迭代路徑呈鋸齒狀。此算法開始時步長較大，當(dāng)其越接近最優(yōu)值點，收斂速度越慢。因此后來提出了（ii），下降算法在文[8]定理2.5.5中給出了收斂性證明。但因每次的迭代方向并未明確，所以下降算法效率不穩(wěn)定。（iii）是文[10]中Bertsekas 和Tsitsiklis 提出的，此算法每次以一個函數(shù)的負梯度方向為迭代方向，這樣對求解大型數(shù)據(jù)集時，理論上增量梯度算法會比普通梯度法有更快速率的收斂。但由于此算法選擇的迭代方向不一定是下降方向，因此不能采用線搜索方法來解決增量梯度算法的問題，而目前主要采用有發(fā)散步長規(guī)則（收斂速度較慢，迭代次數(shù)較大）和恒定步長準則（如果恒定正步長α選擇太大，算法產(chǎn)生的點列{xk}可能不收斂;如果恒定正步長α選擇太小，又會增加迭代的次數(shù)）。結(jié)合上述算法，本文提出了求解問題（1）的一種特殊下降算法——分裂梯度法。這一算法主要思想是，根據(jù)目標函數(shù)中一個函數(shù)的負梯度選擇下降方向dk。具體來說，分裂梯度法主要選擇與上一迭代方向dk成銳角且為下降方向的一個負梯度作為迭代方向;否則，選取-SymbolQC@

f（xk）為當(dāng)前下降方向。所以，分裂梯度法本質(zhì)上就是一種下降算法。只是上文提到的下降算法，其迭代方向并沒有具體給出，而分裂梯度法則針對問題（1）明確給定了每次迭代的下降方向。通過數(shù)值實驗可以看出，與上文提到的三種算法比較，大多數(shù)情況下采用分裂梯度算法更有效，更穩(wěn)健。

1"預(yù)備知識

在本文，我們將用到以下基本知識（見文[8，10]）。

設(shè)函數(shù)f∶Rn→R為一階連續(xù)可微函數(shù)。對任意x∈Rn，f在點x的梯度記為SymbolQC@

f（x）。

定義1"設(shè)DRn是Rn一個子集，若存在常數(shù)Lgt;0，使得

SymbolQC@

f（x1）-SymbolQC@

f（x2）≤Lx1-x2， x1，x2∈D

成立，則稱SymbolQC@

f在D上滿足模為L的Lipschitz連續(xù)。

定義2"設(shè)DRn是Rn一個子集，若對任意給定一個正數(shù)εgt;0，存在一個實數(shù)δgt;0，使得對任意的x1，x2∈D，且滿足x1-x2lt;δ，總有

SymbolQC@

f（x1）-SymbolQC@

f（x2）lt;ε，

則稱SymbolQC@

f在D上一致連續(xù)。

顯然若SymbolQC@

f滿足Lipschitz連續(xù)，則SymbolQC@

f一定滿足一致連續(xù)，反之不然。

2"算法及收斂性定理

考慮無約束優(yōu)化問題

minx∈Rnf（x）∶=∑i∈Ifi（x），（2）

其中I={1，2，…，m}，fi（i∈I）∶Rn→R為一階連續(xù)可微函數(shù)。

我們提出了下述算法。

算法2.1"（分裂梯度法）

Step 0"取x0∈Rn，令k=0，εgt;0，σ∈（0，1），dk=-SymbolQC@

f1（xk）。

Step 1"如果SymbolQC@

f（xk）≤ε，算法終止;否則，轉(zhuǎn)下一步。

Step 2"令i=kmodm+1，計算gk=-SymbolQC@

fi（xk），如果

〈-SymbolQC@

f（xk），gk〉SymbolQC@

f（xk）·gkgt;σ"且"〈dk，gk〉dk·gkgt;σ，（3）

令dk=gk;否則，令dk=-SymbolQC@

f（xk）。

Step 3"選取步長因子αkgt;0。令xk+1=xk+αkdk，k=k+1，轉(zhuǎn)Step 1。

對于上述算法，在Step 3 中可選取的步長{αk}的準則主要有以下幾種：

（1）精確一維搜索準則：αk=argminα≥0f（xk+αdk） 。

（2）Wolfe步長準則：令ω∈0，12，αk滿足：

f（xk）+（1-ω）αkSymbolQC@

f（xk）Tdk≤f（xk+αkdk）

≤f（xk）+ωαkSymbolQC@

f（xk）Tdk。

（3）Goldstein步長準則：αk同時滿足下列條件：

f（xk+αkdk）≤f（xk）+σ1αkSymbolQC@

f（xk）Tdk和

SymbolQC@

f（xk+αkdk）Tdk≥σ2SymbolQC@

f（xk）Tdk，

其中，0lt;σ1lt;σ2lt;1。

（4）Armijo步長準則：αk=βγnk，其中nk為滿足下式的最小非負整數(shù)n：

f（xk+βγndk）≤f（xk）+δβγnSymbolQC@

f（xk）Tdk，

其中βgt;0，δ，γ∈（0，1）為常數(shù)。

上述四種步長準則在下降算法中皆有相似的收斂結(jié)果[8-9]。本文主要以Armijo步長規(guī)則為例，來比較算法2.1與已有算法的收斂效率，包括（隨機）下降算法、最速下降算法、增量梯度法?，F(xiàn)在我們首先給出下降算法和最速下降算法（可參考文[8-9]）。

下降算法。設(shè)當(dāng)前迭代點為xk（k∈N），選擇迭代方向dk滿足dkTSymbolQC@

f（xk）lt;0，利用Armijo步長準則產(chǎn)生步長因子αk。令xk+1=xk+αkdk。

最速下降算法。設(shè)當(dāng)前迭代點xk（k∈N），選擇迭代方向dk=-SymbolQC@

f（xk），利用Armijo步長準則產(chǎn)生步長因子αk。令xk+1=xk+αkdk。

注2.1"（i）由算法2.1中Step 2 的（3）式中的第一個式子可以看出算法2.1本質(zhì)上就是下降算法，算法2.1只是在選取迭代方向時選取單個函數(shù)fi的負梯度或函數(shù)f的負梯度作為下降方向。（ii）算法2.1中Step 2 的（3）式中的第二個式子則是為避免迭代路徑呈鋸齒形，提高收斂速率。

對于下降算法，有以下收斂性結(jié)果，可在[8]定理2.5.5中找到。

命題2.1"假設(shè)目標函數(shù)f在

Rn上一階連續(xù)可微有下界，其梯度函數(shù)SymbolQC@

f在水平集

{xf（x）≤f（x0）}上一致連續(xù)。設(shè){xk}是下降算法采用Armijo步長規(guī)則產(chǎn)生的以x0為初始點的序列。若存在0lt;θ≤π2，使搜索方向dk與-SymbolQC@

f（xk）的夾角θk滿足

θk≤π2-θ，k∈N，（4）

則有 limk→SymboleB@

‖SymbolQC@

f（xk）‖=0 。

下面給出算法2.1（分裂梯度法）的收斂性結(jié)果。

定理2.1"假設(shè)函數(shù)fi（i∈I）在

Rn上一階連續(xù)可微有下界，梯度函數(shù)SymbolQC@

f在水平集{xf（x）≤f（x0）}上一致連續(xù)。設(shè){xk}是算法2.1采用Armijo步長規(guī)則產(chǎn)生的以x0為初始點的序列，則有 limk→SymboleB@

SymbolQC@

f（xk）=0 。

證明"設(shè)k∈N，θk是搜索方向dk與-SymbolQC@

f（xk）的夾角。由Step 2 的（3）式中的第一個式子得

cosθk=〈-SymbolQC@

f（xk），dk〉SymbolQC@

f（xk）·dk

=〈-SymbolQC@

f（xk），gk〉SymbolQC@

f（xk）·gk≥σ〈-SymbolQC@

f（xk），-SymbolQC@

f（xk）〉SymbolQC@

f（xk）·-SymbolQC@

f（xk）=1，

其中σ由算法2.1的Step 0給出，所以cosθk≥σ，θk≤arccosσ。因此（4）式成立，其中θ=π2-arccosσ。所以由命題2.1，有 limk→SymboleB@

SymbolQC@

f（xk）=0。

為求解問題（2），有每次考慮一個函數(shù)的增量梯度算法[10-11]，算法迭代格式如下。

增量梯度算法[4]。設(shè)當(dāng)前迭代點xk（k∈N）。其迭代方式為：xk+1=φm，其中φi=φi-1+αkSymbolQC@

fi（φi-1）（i=1，2，…，m），φ0=xk。其中步長αk滿足

αkgt;0，∑SymboleB@

k=0αk=SymboleB@

，∑SymboleB@

k=0α2klt;SymboleB@

。（5）

關(guān)于增量梯度算法，有以下收斂性結(jié)果（見文[10] Proposition 2）。設(shè)函數(shù)fi（i∈I）在Rn上一階連續(xù)可微有下界，梯度SymbolQC@

fi在Rn上滿足模為L的Lipschitz條件， 并且存在實數(shù)C，Dgt;0，使得SymbolQC@

fi（x）≤C+DSymbolQC@

f（x），（x∈Rn，i∈I）。設(shè){xk}為增量梯度算法產(chǎn)生的點列，則有l(wèi)imk→SymboleB@

SymbolQC@

f（xk）=0。

3"數(shù)值算例

在這一節(jié)，我們將用算法2.1求解幾個具體優(yōu)化問題（2），并通過數(shù)值計算結(jié)果比較該算法和最速下降算法、隨機下降算法以及增量梯度算法的收斂效率（這里的隨機下降算法是指隨機產(chǎn)生下降方向的下降算法）。特別地，我們應(yīng)用算法2.1求解例3.2和例3.3中的穩(wěn)健估計問題和源定位問題，這兩個算例均來自參考文獻[4]，可以轉(zhuǎn)化為問題（2）求解，有重要的實際應(yīng)用背景。在所有算例中算法計算精度均取為ε=10-6，最大迭代次數(shù)設(shè)為5000次，編程軟件為Matlab 2016a[12]。最速下降算法、隨機下降算法和算法2.1均采用Armijo步長規(guī)則，參數(shù)選取為β=1，δ=0.4，γ=0.5;而增量梯度算法采用滿足（5）式的迭代步長αk=1k+k0，（k0∈N） （例3.1選擇的k0較大，是為保證算法最終得到較好的運行結(jié)果;例3.2和例3.3中取k0=1）。

例3.1"設(shè)f1，f2∶R2→R分別定義為：對（x1，x2）∈R2，

f1（x1，x2）=（1-x1）2，f2（x1，x2）=100（x2-x12）2。

顯然函數(shù)fi（i=1，2）一階連續(xù)可微，容易驗證問題minx∈R2f（x）=f1（x）+f2（x）有唯一最優(yōu)解x*=（1，1），最優(yōu)值f（x*）=0。例3.1的實驗結(jié)果見表1。

例3.2"穩(wěn)健估計問題[4]。傳感器網(wǎng)絡(luò)是部署大量低成本傳感器來密集監(jiān)視某個區(qū)域的某種性態(tài)。由于低成本傳感器的可靠性有限，系統(tǒng)必須設(shè)計成對單個傳感器的可能故障具有魯棒性。這意味著在評估任務(wù)中一些傳感器會產(chǎn)生不合理的測量結(jié)果，即異常值。在文獻[13]中，作者建議使用穩(wěn)健

的統(tǒng)計數(shù)據(jù)來減輕數(shù)據(jù)中的異常值的影響（見[14]或[15]）。穩(wěn)健估計使用的是賦予異常值較少權(quán)重的目標函數(shù)，可以實現(xiàn)此目的的一個常用函數(shù)是“Fair”函數(shù)g∶R→R[16]，定義為

g（x）=c2xc-ln1+xc，x∈R。（6）

與參考文獻[4]類似，我們模擬了一個用于測量污染水平的傳感器網(wǎng)絡(luò)，并假設(shè)一定比例的傳感器損壞且提供不合理的測量結(jié)果。每個傳感器都收集污染水平的單個噪聲測量值，通過最小化目標

函數(shù)（7）來確定平均污染水平的估計值。這里研究的穩(wěn)健估計問題轉(zhuǎn)化為優(yōu)化問題

minx∈R

fx=∑mi=1fi（x），（7）

其中fi∶R→R定義為

fi（x）=1mg（x-yi）， x∈R。

其中yi是由傳感器i采集的測量值。容易驗證函數(shù)fi（i=1，2，…，m）一階連續(xù)可微?，F(xiàn)假設(shè)本次模擬有20個傳感器，即m=20。為了反映傳感器故障的可能性，一半的樣本按均值為10，方差為1的高斯分布生成，另一半是按均值為10，方差為10高斯分布生成。其中（6）式中的系數(shù)c=10。例3.2的實驗結(jié)果見表2，其中增量梯度算法的步長αk=1k+1。

例3.3"源定位問題[4]。設(shè)傳感器i等距分布在100*100場的空間位置ri∈R2上，每個傳感器收集從源點x*∈R2所發(fā)射聲信號的噪聲測量yi。源定位問題就是根據(jù)收集的信號{yi}，找到源的位置x*。基于遠場假設(shè)和各向同性聲波傳播模型[13，17-18]，源位置估計問題可歸結(jié)為非線性最小二乘問題：

minx∈R2f（x）=∑mi=1fi（x），

其中fi∶R2→Ri=1，2，…，m定義為

fi（x）=（yi-g（ri-x2））2， x∈R2。（8）

函數(shù)g：R→R定義為

g（z）=Az，

z≥Aε

2ε-zε21000，zlt;Aε，z∈R。 （9）

在上式中，A和ε是表示與信號強度相關(guān)的常數(shù)。容易驗證函數(shù)fi（i=1，2，…，m）一階連續(xù)可微。在我們的數(shù)值模擬實驗中，源點x*=（60，60），A=1000[4]，ε=1，m=16。取ri為100*100網(wǎng)格上均勻分布的16個點，yi是根據(jù)高斯分布產(chǎn)生的，其均值為g（ri-x*2），方差為1。例3.3的實驗結(jié)果見表3，其中增量梯度算法中的步長αk=1k+1。

注3.1"（i）算法2.1與最速下降算法比較。由例3.1和例3.3的數(shù)值結(jié)果可知，取不同的初始點，算法2.1的收斂效率明顯都高于最速下降算法。在例3.2中，算法2.1的迭代次數(shù)大于最速下降算法的迭代次數(shù)，分析其原因為自變量x是一維的，采用最速下降算法時，迭代路徑不會出現(xiàn)鋸齒形，數(shù)值結(jié)果表明最速下降算法的收斂速率快于其他算法。

（ii）算法2.1與下降算法比較。由例3.1～3.3的數(shù)值結(jié)果可知，對于不同的初始點，算法2.1的收斂效率平均優(yōu)于隨機下降算法，并且可以看出隨機下降算法的迭代次數(shù)并不穩(wěn)定。

（iii）算法2.1與增量梯度算法比較。對比同樣考慮由多個函數(shù)之和構(gòu)成的增量梯度算法，由例3.1～3.3的數(shù)值結(jié)果表明在相同精度要求下，算法2.1的迭代次數(shù)明顯少于增量梯度算法，且在控制迭代步數(shù)（5000步）內(nèi)都達到計算精度要求，而增量梯度算法均運行5000次，未達到計算精度要求。

參考文獻：

[1]錢偉懿，張洵. 一種改進的粒子群優(yōu)化算法[J]. 渤海大學(xué)學(xué)報（自然科學(xué)版），2017，38（2）：97-103.

[2]Solodov M. V. Incremental gradient algorithms with step sizes bounded away from zero[J]. Comput. Optim. Appl.， 1998， 11（1）：23-35.

[3]Grippo L. A class of unconstrained minimization methods for neural network training[J]. Optimization Methods and Software，1994， 4（2）：135-150.

[4]Blatt D， Hero A，Gauchman H. A convergent incremental gradient method with a constant step size[J]. SIAM Journal on Optimization，2007，18（1）：29-51.

[5]Bertsekas D P. A new class of incremental gradient methods for least squares problems[J]. SIAM Journal on Optimization， 1997， 7（4）：913-926.

[6]Tseng P. "An incremental gradient（-projection） method with momentum term and adaptive step size rule[J]. SIAM Journal on Optimization， 1998，8：506-531.

[7]Defazio A， Bach F， Lacoste￣Julien S. A fast incremental gradient method with support for non￣strongly convex composite objectives[J]. International Conference on Neural Information Processing Systems， 2014，1：1646-1654.

[8]袁亞湘，孫文瑜. 最優(yōu)化理論與方法[M].北京：科學(xué)出版社， 2007.

[9]王宜舉，修乃華. 非線性最優(yōu)化理論與方法[M].北京：科學(xué)出版社， 2015.

[10]Bertsekas D P，Tsitsiklis J N. Gradient convergence in gradient methods with errors[J]. SIAM Journal on Optimization，2000， 10（3）：627-642.

[11]Bertsekas D P. 非線性規(guī)劃[M]. 2版.北京：清華大學(xué)出版社， 2013.

[12]Trefethen. Spectral Methods in MATLAB[M]. New York： SIAM， 2000.

[13]Rabbat M G，Nowak R D. Decentralized Source Localization and Tracking[C]. Montreal： IEEE International Conference on Acoustics， 2004.

[14]Huber P. Robust Statistics[M]. New York： John Wiley amp; Sons， 1981.

[15]Polyak B T. Introduction to Optimization[M]. New York： Optimization Software Income，1987.

[16]Rey W J J. Introduction to Robust and Quasi￣robust Statistical Methods[M]. Berlin： Springer Verlag， 1983.

[17]Sheng X H， Hu Y H. Information Processing in Sensor Networks[M]. California： Springer， 2003.

[18]Chen J C， Yao K， Hudson R E. Source localization and beam forming[J]. IEEE Signal Processing Magazine， 2002（19）：30-39.

（責(zé)任編輯：曾"晶）

A Special Descent Algorithm——Split Gradient Method

QIAN Xiaohui，WANG Xiangmei*

（College of Mathematics and Statistics，Guizhou University，Guiyang 550025，China）

Abstract：

The common methods to solve unconstrained optimization problems include the descent algorithm， Newton method and the conjugate gradient method， etc. When the objective function is the sum of several smooth functions， some authors propose and study the incremental gradient algorithm. The algorithm cyclically select the negative gradient of a single function as the iteration direction， which is not necessarily a descent direction. Therefore， incremental gradient algorithm is not a descent algorithm in general. In this paper， the split gradient method was proposed which combines the ideas of the descent algorithm and the incremental gradient algorithm. Finally， some numerical experiments "were provided to compare the split gradient method with the steepest descent algorithm， the random descent algorithm and the incremental gradient algorithm， respectively. The numerical results show that the split gradient method is more effective than the others for some optimization problems.

Key words：

unconstrained optimization; descent algorithm; incremental gradient method; split gradient method; Armijo step rule