梯度的方向與等值面垂直,並且指向函數值提升的方向。
二次收斂是指一個算法用於具有正定二次型函數時,在有限步可達到它的極小點。二次收斂與二階收斂沒有盡然聯系,更不是一回事,二次收斂往往具有超線性以上的收斂性。一階收斂不一定是線性收斂。
解釋一下什么叫正定二次型函數:
n階實對稱矩陣Q,對於任意的非0向量X,如果有XTQX>0,則稱Q是正定矩陣。
對稱矩陣Q為正定的充要條件是:Q的特征值全為正。
二次函數,若Q是正定的,則稱f(X)為正定二次函數。
黃金分割法
黃金分割法適用於任何單峰函數求極小值問題。
求函數在[a,b]上的極小點,我們在[a,b]內取兩點c,d,使得a<c<d<b。並且有
1)如果f(c)<f(d),則最小點出現在[a,d]上,因此[a,d]成為下一次的搜索區間。
2)如果f(c)>f(d),則[c,b]成為下一次的搜索區間。
假如確定了[a,d]是新的搜索區間,我們並不希望在[a,d]上重新找兩個新的點使之滿足(1)式,而是利用已經抗找到有c點,再找一個e點,使滿足:
可以解得r=0.382,而黃金分割點是0.618。
練習:求函數f(x)=x*x-10*x+36在[1,10]上的極小值。
#include<stdio.h> #include<math.h> #include<limits.h> double func(double x){ return x*x-10*x+36; } void main(){ double zeta=0.001; double a=1.0,b=10.0; double t1=a-1; double t2=b+1; double v1=0.0; double v2=0.0; double min_value=INT_MAX; int iteration=0; while(++iteration){ if(t1==a-1){ t1=a+0.382*(b-a); v1=func(t1); } if(t2==b+1){ t2=a+0.618*(b-a); v2=func(t2); } if(v1<v2){ min_value=v1; b=t2; t2=t1; v2=v1; t1=a-1; } else{ min_value=v2; a=t1; t1=t2; v1=v2; t2=b+1; } if(fabs(b-a)<zeta) break; printf("當前極小值%f\n",min_value); } printf("迭代次數%d\n",iteration); printf("極小值%f\n",min_value); }
最速下降法
泰勒級數告訴我們:
其中Δx可正可負,但必須充分接近於0。
X沿D方向移動步長a后,變為X+aD。由泰勒展開式:
目標函數:
a確定的情況下即最小化:
向量的內積何時最小?當然是兩向量方向相反時。所以X移動的方向應該和梯度的方向相反。
接下來的問題是步長a應該怎么定才能使迭代的次數最少?
若f(X)具有二階連續偏導,由泰勒展開式可得:
H是f(X)的Hesse矩陣。
可得最優步長:
g是f(X)的梯度矩陣。
此時:
可見最速下降法中最優步長不僅與梯度有關,而且與Hesse矩陣有關。
練習:求函數f(x1,x2)=x1*x1+4*x2*x2在極小點,以初始點X0=(1,1)T。
#include"matrix.h" #include<iostream> #include<iomanip> #include<cmath> #include<limits> #include<cassert> using namespace std; const int SIZE=2; const double ZETA=0.001; inline double func(Matrix<double> &X){ assert(SIZE==X.getRows()); double x1=X.get(0,0); double x2=X.get(1,0); return x1*x1+4*x2*x2; } inline Matrix<double> gradient(Matrix<double> &X){ assert(SIZE==X.getRows()); double x1=X.get(0,0); double x2=X.get(1,0); Matrix<double> rect(SIZE,1); rect.put(0,0,2*x1); rect.put(1,0,8*x2); return rect; } inline Matrix<double> Hesse(Matrix<double> &X){ Matrix<double> rect(SIZE,SIZE); rect.put(0,0,2); rect.put(0,1,0); rect.put(1,0,0); rect.put(1,1,8); return rect; } int main(int agrc,char *argv[]){ Matrix<double> X(SIZE,1); X.put(0,0,1); X.put(1,0,1); int iteration=0; double value=func(X); double newValue=numeric_limits<double>::max(); while(++iteration){ Matrix<double> G=gradient(X); double factor=((G.getTranspose()*G).get(0,0))/((G.getTranspose()*Hesse(X)*G).get(0,0)); for(int i=0;i<G.getRows();++i) for(int j=0;j<G.getColumns();++j) G.put(i,j,G.get(i,j)*factor); Matrix<double> newX=X-G; //cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl; newValue=func(newX); if(fabs(newValue-value)<ZETA) break; else{ X=newX; value=newValue; } //cout<<"本次迭代找到的極小值"<<value<<endl; } cout<<"迭代次數"<<iteration<<endl; cout<<"極小值"<<value<<endl; return 0; }
梯度下降法開始的幾步搜索,目標函數下降較快,但接近極值點時,收斂速度就比較慢了,特別是當橢圓比較扁平時,收斂速度就更慢了。
另外最速下降法是以函數的一次近似提出的,如果要考慮二次近似,就有牛頓迭代法。
牛頓迭代法
在點Xk處對目標函數按Taylar展開:
令
得
即
可見X的搜索方向是,函數值要在此方向上下降,就需要它與梯度的方向相反,即
。所以要求在每一個迭代點上Hesse矩陣必須是正定的。
練習:求的極小點,初始點取X=(0,3)。
#include"matrix.h" #include<iostream> #include<iomanip> #include<cmath> #include<limits> #include<cassert> using namespace std; const int SIZE=2; const double ZETA=0.001; inline double func(Matrix<double> &X){ assert(SIZE==X.getRows()); double x1=X.get(0,0); double x2=X.get(1,0); return pow(x1-2,4)+pow(x1-2*x2,2); } inline Matrix<double> gradient(Matrix<double> &X){ assert(SIZE==X.getRows()); double x1=X.get(0,0); double x2=X.get(1,0); Matrix<double> rect(SIZE,1); rect.put(0,0,4*pow(x1-2,3)+2*(x1-2*x2)); rect.put(1,0,-4*(x1-2*x2)); return rect; } inline Matrix<double> Hesse(Matrix<double> &X){ Matrix<double> rect(SIZE,SIZE); double x1=X.get(0,0); double x2=X.get(1,0); rect.put(0,0,12*pow(x1-2,2)+2); rect.put(0,1,-4); rect.put(1,0,-4); rect.put(1,1,8); return rect; } int main(int agrc,char *argv[]){ Matrix<double> X(SIZE,1); X.put(0,0,0); X.put(1,0,3); int iteration=0; double value=func(X); double newValue=numeric_limits<double>::max(); while(++iteration){ Matrix<double> G=gradient(X); Matrix<double> H=Hesse(X); Matrix<double> newX=X-H.getInverse()*G; cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl; newValue=func(newX); if(fabs(newValue-value)<ZETA) break; else{ X=newX; value=newValue; } cout<<"本次迭代找到的極小值"<<value<<endl; } cout<<"迭代次數"<<iteration<<endl; cout<<"極小值"<<value<<endl; return 0; }
牛頓法是二次收斂的,並且收斂階數是2。一般目標函數在最優點附近呈現為二次函數,於是可以想像最優點附近用牛頓迭代法收斂是比較快的。而在開始搜索的幾步,我們用梯度下降法收斂是比較快的。將兩個方法融合起來可以達到滿意的效果。
收斂快是牛頓迭代法最大的優點,但也有致命的缺點:Hesse矩陣及其逆的求解計算量大,更何況在某個迭代點Xk處Hesse矩陣的逆可能根本就不存在(即Hesse矩陣奇異),這樣無法求得Xk+1。
擬牛頓法
Hesse矩陣在擬牛頓法中是不計算的,擬牛頓法是構造與Hesse矩陣相似的正定矩陣,這個構造方法,使用了目標函數的梯度(一階導數)信息和兩個點的“位移”(Xk-Xk-1)來實現。有人會說,是不是用Hesse矩陣的近似矩陣來代替Hesse矩陣,會導致求解效果變差呢?事實上,效果反而通常會變好。
擬牛頓法與牛頓法的迭代過程一樣,僅僅是各個Hesse矩陣的求解方法不一樣。
在遠離極小值點處,Hesse矩陣一般不能保證正定,使得目標函數值不降反升。而擬牛頓法可以使目標函數值沿下降方向走下去,並且到了最后,在極小值點附近,可使構造出來的矩陣與Hesse矩陣“很像”了,這樣,擬牛頓法也會具有牛頓法的二階收斂性。
對目標函數f(X)做二階泰勒展開:
兩邊對X求導
當X=Xi時,有
這里我們用Hi來代表在點Xi處的Hesse矩陣的逆,則
(5)式就是擬牛頓方程。
下面給出擬牛頓法中的一種--DFP法。
令
我們希望Hi+1在Hi的基礎上加一個修正來得到:
給定Ei的一種形式:
m和n均為實數,v和w均為N維向量。
(6)(7)聯合起來代入(5)可得:
下面再給一種擬牛頓法--BFGS算法。
(8)式中黑色的部分就是DFP算法,紅色部分是BFGS比DFP多出來的部分。
BFGS算法不僅具有二次收斂性,而且只有初始矩陣對稱正定,則BFGS修正公式所產生的矩陣Hk也是對稱正定的,且Hk不易變為奇異,因此BFGS比DFP具有更好的數值穩定性。
共軛方向法
最速下降法有鋸齒現像,收斂速度慢;而牛頓法需要計算Hesse矩陣而計算量大。共軛方向法收斂速度界於兩者之間,具有二次收斂性。共軛方向法屬於效果好而又實用的方法。
由於一般目標函數在最優點附近呈現為二次函數,因此可以設想一個算法對於二次函數比較有效,就可能對一般函數也有較好效果。共軛方向法是在研究對稱正定二次函數的基礎上提出來的。
則稱兩個向量P0和P1為Q的共軛向量。當Q為單位向量時,有,所以“共軛”是“正交”的推廣。
對於二次正定函數,從任意點X0出發,沿任意下降方向P0作直線搜索得到X1,再從X1出發,沿與P0共軛的方向P1作直線搜索,即可得到f(X)的極小點。
當一組向量Pi(i=1,2,...,n-1)為Q共軛時,從任意點出發,依次沿P0,P1,P2,...,Pn-1方向作下述算法的直線搜索,經過n次迭代必定收斂於正定二次函數的極小點。
為確定最優步長tk,令
現在問題是如何產生一組關於Q共軛的向量?這里一種叫作Gram-Schmidt的方法。
取線性無關的向量組V0,V1,...,Vn-1,例如取n個坐標軸的單位向量。
取P0=V0.
上面的方法都是針對目標函數為正定二次函數的,對於一般非二次函數,可以通過二次近似。
這就是f(X)在極小點X*處的近似,是Hesse矩陣,相當於Q,由於X*未知,但當X0充分接近於X*時,可用
近似代替
,從而構造共軛向量。
理論與實踐證明,將二次收斂算法用於非二次的目標函數,亦有很好的效果,但迭代次數不一定保證有限次,即對非二次n維目標函數經n步共軛方向一維搜索不一定就能達到極小點。在這種情況下,為了找到極小點,可用泰勒級數將該函數在極小點附近展開,略去高於二次的項之后即可得該函數的二次近似。實際上很多的函數都可以用二次函數很好地近似,甚至在離極小點不是很近的點也是這樣。故用二次函數近似代替非二次函數來處理的方法不僅在理論分析上是重要的,而且在工程實際應用中也是可取的。
共軛梯度法
共軛梯度法是共軛方向法的一種延伸,初始共軛向量P0由初始迭代點X0處的負梯度-g0來給出。以后的Pk由當前迭代點的負梯度與上一個共軛向量的線性組合來確定:
對於非二次函數的優化問題,迭代次數不止n次,但共軛方向只有n個。當迭代n次后,可以把Pn重新置為最開始的P0,其他的變量按原方法更新。