## 無約束最優化方法 最優化(Optimization)是工程技術、經濟管理、 科學研究、社會生活中經常遇到的問題, 如：結構設計，資源分配，生產計劃，運輸方案等等。 最優化: 在一定條件下，尋求使目標最大(小)的決策解決優化問題的手段 ? 經驗積累，主觀判斷 ? 作試驗，比優劣 ? 建立數學模型，求解最優策略 ### 優化問題的數學模型 無約束優化問題標準形式： 求n維決策變量x=[x1,x2,…,xn]T在可行域Ω中 使得目標函數f(x)→min,x?Ω?Rn ?  **可行解**只滿足(2)，**最優解**滿足(1)(2) **無約束優化**只有(1)，**約束優化**需要(1)(2) ### 無約束優化的基本方法 ### 求解無約束優化的基本思路：[下山法](http://blog.csdn.net/u013007900/article/details/43525733) #### 基本思想： 和爬山法基本思想一樣，在Rn中某一點，確定一個**搜索放向**然后沿著該方向的**移動步長**，得到使目標函數下降的最新點。 #### 迭代步驟： Step 1 初始化：初始點x0，終止準則等 Step 2 迭代改進：方向dk，步長αk xk+1=xk+(αd)k?且?f(xk+1)<f(xk) Step 3 終止檢驗：得到近優解或k+1?k轉2> 選擇dk,αk使得f下降更快?利用不同的算法 ### 最速下降法（梯度法） 泰勒級數：  其中Δx可正可負，但必須充分接近于0。 X沿D方向移動步長α后，變為X+αD。由泰勒展開式： 目標函數： α確定的情況下即最小化： 向量的內積何時最小？當然是兩向量方向相反時。所以X移動的方向應該和梯度的方向相反。 **下降方向：**?fT(xk)dk<0 **最速下降方向：**dk=??f(xk)負梯度方向 **迭代改進形式：**xk+1=xk??f(xk) 算法特點：在初始狀態改進較快，在最優解附近改進緩慢。 求函數f(x1,x2)=x21+4x22在極小點，以初始點X0=(1,1)T。 ~~~ #include"matrix.h"//大家自己YY一個矩陣類吧 #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; } ~~~ > 這份算法使用的是利用Hesse矩陣和泰勒公式計算最優步長 若f(X)具有二階連續偏導，由泰勒展開式可得：  H是f(X)的Hesse矩陣。  可得最優步長：  g是f(X)的梯度矩陣。 此時：  可見最速下降法中最優步長不僅與梯度有關，而且與Hesse矩陣有關。 ### 牛頓迭代法 將f(xk+1)在xk點作泰勒展開至二階項，用d替代dk，有 f(xk+1)=f(xk+d)=f(xk)+?fT(xk)d+12dT?2f(xk)d 求d使得f(xk+1)極小?右端對d的導數為0??f(xk)+?2f(xk)d=0 可得到> **牛頓方程**?2f(xk)d?=???f(xk) **牛頓方向**dk?=??(?2f(xk))?1?f(xk) **迭代格式**xk+1?=?xk?(?2f(xk))?1?f(xk) 其中?2f(xk)為[Hessen矩陣](http://baike.baidu.com/link?url=E2wdH1BhTKFos1JYd_NhHvkmGAyABsHJDIgSyygvzTuoAs41tvb4OGzzWThbUMwobdBAEPpS7nzia-GQUHHf7a) 可見x的搜索方向是?(?2f(xk))?1?f(xk)，函數值在此方向上下降，就需要它與梯度方向相反，即??f(xk)(?2f(xk))?1?f(xk)<0，所以要求在每一個迭代點上Hessian矩陣必須是正定的。 算法特點： 牛頓法是二次收斂的，并且收斂階數是2。一般目標函數在最優點附近呈現為二次函數，于是可以想像最優點附近用牛頓迭代法收斂是比較快的。而在開始搜索的幾步，我們用梯度下降法收斂是比較快的。將兩個方法融合起來可以達到滿意的效果。 收斂快是牛頓迭代法最大的優點，但也有致命的缺點：Hessian矩陣及其逆的求解計算量大，更何況在某個迭代點Xk處Hessian矩陣的逆可能根本就不存在（即Hessian 矩陣奇異），這樣無法求得Xk+1。 求f(x1,x2)=(x1?2)4+(x1?2x2)2的極小點，初始點取X=(0,3)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 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; } ~~~ ### 擬牛頓法 Hessian矩陣在擬牛頓法中是不計算的，擬牛頓法是構造與Hessian矩陣相似的正定矩陣，這個構造方法，使用了目標函數的梯度（一階導數）信息和兩個點的“位移”（Xk?Xk?1）來實現。 有人會說，是不是用Hessian矩陣的近似矩陣來代替Hessian矩陣，會導致求解效果變差呢？事實上，效果反而通常會變好。 擬牛頓法與牛頓法的迭代過程一樣，僅僅是各個Hessian矩陣的求解方法不一樣。 在遠離極小值點處，Hessian矩陣一般不能保證正定，使得目標函數值不降反升。而擬牛頓法可以使目標函數值沿下降方向走下去，并且到了最后，在極小值點附近，可使構造出來的矩陣與Hessian矩陣“很像”了，這樣，擬牛頓法也會具有牛頓法的二階收斂性。 對目標函數f(X)做二階泰勒展開： f(x)=f(xk+1)+(x?xk+1)T?f(xk+1)+12(x?xk+1)TG(xk+1)(x?xk+1) 其中?2f(xk)??G(xk+1) 上式兩邊對x求導： ?f(x)=?f(xk+1)+G(xk+1)(x?xk+1) 當x=xk的時候，有 G?1(xk+1)[?f(xk+1)??f(xk)]=xk+1?xk 這里我們用Hk表示在點xk處的Hessian矩陣的**逆矩陣** Hk+1[?fd(xk+1)??f(xk)]=xk+1?xk?......(5) (5)式就是**擬牛頓方程** 下面提供兩種擬牛頓法  ### 非線性最小二乘(Least Square)擬合 問題： 給定(ti,yi),i=1,…n, 擬合一個函數y=f(t,x)， 其中x為待定的參數向量, f 對x非線性。 記誤差ri(x)=yi?f(ti,x)?????r(x)=(r1(x),…,rn(x))T 優化模型： minxR(x)=rT(x)r(x) 記r(x)的雅克比行列式為J(x)=(?ri?xj)n×m ?R=2J(x)Tr(x)?2R=2J(x)TJ(x)+2SS=∑ni=1ri(x)?2ri(x)?2ri(x)=(?2ri?xk?xl)n×m 牛頓法要計算Hessian矩陣，其中S計算量大 若f 對x線性, 則化為線性最小二乘擬合, 此時 S=0 特定算法考慮如何忽略或近似矩陣S #### Gauss-Newton算法(GN)：忽略矩陣S ?R=2J(x)Tr(x)?2R=2J(x)TJ(x) 牛頓方程：?2f(xk)d=??f(xk) f用R代替，下降方向dk滿足 J(x)TJ(x)dk=?J(x)Tr(x) GN算法收斂性特點: S=∑ni=1ri(x) 收斂性依賴f對x的線性程度，即偏差r的大小（”小殘差問題“） #### LM算法: Levenbery (1944), Marquardt (1963) G-N算法修正 J(x)TJ(x)dk=?J(x)Tr(x) ? (J(x)TJ(x)+αkI)dk=?J(x)Tr(x) 其中αk>0為修正參數，防止JTJ出現病態。 LM算法特點 dk位于GN方向（αk很小）和負梯度方向(αk很大)之間 LM算法的另一種觀點 min∥r(xk)+J(xk)(x?xk)∥2 s.t.?????∥x?xk∥2≤hk????(?)置信域 其解為 (J(x)TJ(x)+αkI)dk=?J(x)Tr(x) 若GN方向滿足約束(*)，αk=0,否則ak>0 置信域算法特點： 先確定最長步長hk，再確定方向dk