錢曉慧 王湘美

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

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

中圖分類號：O224

文獻標識碼： A

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

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

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

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

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

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

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

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

1?預備知識

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

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

f（x）。

定義1?設DRn是Rn一個子集，若存在常數L>0，使得

SymbolQC@

f（x1）-SymbolQC@

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

成立，則稱SymbolQC@

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

定義2?設DRn是Rn一個子集，若對任意給定一個正數ε>0，存在一個實數δ>0，使得對任意的x1，x2∈D，且滿足x1-x2<δ，總有

SymbolQC@

f（x1）-SymbolQC@

f（x2）<ε，

則稱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ù)可微函數。

我們提出了下述算法。

算法2.1?（分裂梯度法）

Step 0?取x0∈Rn，令k=0，ε>0，σ∈（0，1），dk=-SymbolQC@

f1（xk）。

Step 1?如果SymbolQC@

f（xk）≤ε，算法終止;否則，轉下一步。

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

fi（xk），如果

〈-SymbolQC@

f（xk），gk〉SymbolQC@

f（xk）·gk>σ?且?〈dk，gk〉dk·gk>σ，（3）

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

f（xk）。

Step 3?選取步長因子αk>0。令xk+1=xk+αkdk，k=k+1，轉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，

其中，0<σ1<σ2<1。

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

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

f（xk）Tdk，

其中β>0，δ，γ∈（0，1）為常數。

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

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

f（xk）<0，利用Armijo步長準則產生步長因子αk。令xk+1=xk+αkdk。

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

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

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

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

命題2.1?假設目標函數f在

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

f在水平集

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

f（xk）的夾角θk滿足

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

則有 limk→SymboleB@

‖SymbolQC@

f（xk）‖=0 。

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

定理2.1?假設函數fi（i∈I）在

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

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

SymbolQC@

f（xk）=0 。

證明?設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），有每次考慮一個函數的增量梯度算法[10-11]，算法迭代格式如下。

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

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

αk>0，∑SymboleB@

k=0αk=SymboleB@

，∑SymboleB@

k=0α2k<SymboleB@

。（5）

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

fi在Rn上滿足模為L的Lipschitz條件， 并且存在實數C，D>0，使得SymbolQC@

fi（x）≤C+DSymbolQC@

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

SymbolQC@

f（xk）=0。

3?數值算例

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

例3.1&nbsp;設f1，f2∶R2→R分別定義為：對（x1，x2）∈R2，

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

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

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

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

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

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

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

minx∈R

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

其中fi∶R→R定義為

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

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

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

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

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

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

函數g：R→R定義為

g（z）=Az，

z≥Aε

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

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

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

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

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

參考文獻：

[1]錢偉懿，張洵. 一種改進的粒子群優(yōu)化算法[J]. 渤海大學學報（自然科學版），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].北京：科學出版社， 2007.

[9]王宜舉，修乃華. 非線性最優(yōu)化理論與方法[M].北京：科學出版社， 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版.北京：清華大學出版社， 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 & 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.

（責任編輯：曾?晶）

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