FR共轭梯度法与拟牛顿法计算机实现及仿真.docx

FR共轭梯度法与拟牛顿法计算机实现及仿真信计 0701 作者：翟铮 学号：5摘要：最速下降法和牛顿法是最基本的无约束最优化方法，但由于最速下降法收敛速度慢而牛顿法虽然收敛速度比 较快，但需要计算目标函数的Hesse矩阵及其逆矩阵，计算量比较大我们这里采用采用一种无需计算二阶导数， 但收敛速度比较快的一种算法，即FR共轭梯度关键词:FR共轭梯度法、拟牛顿法、黄金分割法、最优一维线性搜索、收敛性一、FR共轭梯度法与拟牛顿法1、共轭梯度法在共轭梯度方向法中，如果取初始的搜索方向d二呵(x )，而以以下各共轭方向d由第k次迭代点的负梯0 0 k度-Yf (x )与已经得到的共轭方向d 的线性组合来确定，这样就构造了一种具体的共轭方向法因为每一个共 k k-1轭方向都依赖于迭代点处的负梯度，所以称之为共轭梯度法(conjugate gradient method)给定初始点x e Rn，取初始搜索方向d =-Vf (x )，从x出发，沿d进行一维搜索，得到迭代点x 1，一下0 0 0 0 0 1按d = -Vf (x )+a d , k = 0,1, — n 一 2k +1 k+1 k k来构造搜索方向，u的选取应该使得d 和d是共轭的。

因为k k + 1 kdT Qd =-Vf (x )TQd +u dTQdk +1 k k +1 k k k k且要使d 和d是共轭的，应有dT Qd二0，所以k + 1 k k +1 k,k = 0,1 …，n — 2Vf (x)TQddTQdkk于是得到n个搜索方向do,耳，…，dn-1如下:Fletcher_Reeves 公式为 ak||Vf (x )||2 k+1 II Vf (x )l|2kfd =-Vf (x )0 0d = -Vf (x ) + a dk +1 k+1 k kVf (xk丿创dTQdkk2、拟牛顿法：牛顿法成功的关键是利用了 Hesse矩阵提供的曲率信息，但计算Hesse矩阵工作量大，并且有的目标函数的 Hesse 矩阵很难计算，甚至不好求出针对这一问题，拟牛顿法比牛顿法更为有效这类算法仅利用目标函数值和一阶导数的信息，构造出目标函数的曲率近似，使方法具有类似牛顿法的收敛速度快的优点一般的拟牛顿法算法为：给出初始点 xo Rn，Bo (或Ho) w E / > 0,k 二 0.如果kjl-e，停止解B d = -g 得搜索方向d （或计算d =-H g .）k k k k k kaa由线性搜索得步长因子，并令x二x d .k k +1 k k k校正B产生B （或校正H产生H ），使得拟牛顿条件成立。

k k +1 k k +1本次实验采用 DFP 校正，校正公式为：s sT H y yTH— + —k—k— k—k—k kk+1 k sT y yT H yk k k k k二、计算机实现步骤由于该问题设计的问题比较复杂，所用变量多，为防止变量之间的冲突，在下面的程序中全部用函数来实现函数一：计算函数在任意点处的梯度值函数名：gradient_my(f,y,num),参数：f:待求梯度函数 y：待求梯度点 num:函数变量数 函数二：黄金分割法求最优步长函数名： golden_search(f,a,b),参数：f：待求梯度函数 a:搜索区间左端点 b:搜索区间右端点函数三：共轭梯度法函数名： conjugate_gradient(f,x0,error),参数：f：待求梯度函数 x0:初始点 error:允许误差函数四：拟牛顿法函数名： quasi_Newton(f,x0,error),参数：f：待求梯度函数 x0:初始点 error:允许误差三、实验结果及仿真第一题、Rosenbrock函数f (x) = 100(x 一x2)2 + (1 一x )2,2 1 1x 二[一1・2,1]t, x* 二[1,1]T, f (x*)二 00设定黄金分割法的步长区间为[一0.1,0.1],搜索精度为 0.001,设定共轭梯度法误差为： 0.001调用函数： conjugate_gradient(f,[-1.2,1],0.001);共 轭梯度 法经过 70次迭代 收 敛效果 图结论：从图上可以看到利用共轭梯度法求解过程中，从初始点向收敛点迭代的步长是阶段性的变化的，期间有的阶段步长很小，阶段间的步长比较大。

此方法在搜索在收敛中比较快，有效的避免了最速下降法的锯齿效应第三题：Wood函数f (x) = 100(x2 - x )2 + (x - 1)2 + (x -1)2 + 90(X2 - X)2 + ...1 2 1 3 3 4+10.1[(x -1)2 + (x -1)2] +19.8(x -1)(x -1)2 4 2 4x —[-3,-1,-3,-1]t ,x * =[1丄1,1]t ,f (x*) = 00解：我们用上述的共轭梯度法求解，导致搜索精度低，达不到函数停止迭代的下限，当搜索次数增大时，导致了迭 代不收敛，因此我们试图改用拟牛顿法来求解，观察拟牛顿法在处理高维函数时的优异性我们用拟牛顿法设定搜索步长区间为[-15,15]，误差精度为：0.001 得到的结果如下：迭代次数XYZW1-3.00000-1.00000-3.00000-1.000002-2.19396-0.86038-2.27451-0.873813-2.44735-0.57344-2.05230-0.547474-3.089655.91018-2.949035.3130851.830036.503150.985705.777276-2.588367.31350-2.552965.748667-2.077204.07977-2.885778.548658-2.080723.93519-2.932658.636719-1.427401.08110-3.2210010.4600910-0.805270.27410-3.062168.88119761.073141.143600.913940.81642771.061641.124780.927290.84102781.044451.096630.947330.87837791.019631.042880.966660.91740801.016301.038850.967760.92561811.010661.031290.970720.94075820.977680.960091.011251.02000830.968220.939311.024781.04712840.988390.978461.009751.02145851.000251.000871.000221.00027861.000201.000420.999780.99957870.999990.999981.000001.00001收敛性与迭代次数的关系拟牛顿法 Wood 函数迭代 87次收敛效果示意图2拟牛 顿法 迭 代87次收 敛 效果 示意 图1210收敛点（1-1-2-3初始占t\\\-410-5第四题：Powell奇异函数:/ (X) = (X + 10X )2 + 5(X 一 X )2 + (X 一 2X )4 + 10(x 一 x )41 2 3 4 2 3 1 4x = [-3,-1,0,1]t , x * = [0,0,0,0]t , / (x *) = 0 0解：设定黄金分割法的步长区间为： [-0.1,0.1]，搜索精度为： 0.0001，共轭梯度法的终止点误差范围为IIV/(x*)11 < 0.0001。

迭代终止点为：[0.0673 -0.0067 0.0332 0.0334]2调用函数： conjugate_gradient(f,[3,-1,0,1],0.0001); 可得结果为：由于函数为四维，我们以下下分别用三个图形来表示收敛效果，分别列出了图 1、X-Y-Z 图 2、X-Y-W 图三 Y-Z-W梯度平方和误差精度为 0.0001时收敛效果收敛点251-0.50.70.60.5初始点2.53图12一1.5一2.5-0.5-3YX0：0-1 0图2共轭梯度法梯度平方和精度为 0.0001 时收敛效果Z图3结论：从上图可以看到，共轭梯度法搜索过程中，初始步长较大，进行搜索，随着距离收敛点的距离面小，步长逐步 变小，向最优解逼近共轭梯度法主程序：function A=conjugate_gradient(f,x0,error)[a,b]=size(x0);initial_gradient=gradient_my(f,x0,b);norm=0;norm0=0;syms step_zzh;A=[];for i=1:bnormO=normO+(initial_gradient(i))入2;endsearch_direction=-initial_gradient;x=x0+step_zzh*search_direction;f_step=subs(f,findsym(f),x);best_step=golden_search(f_step,-0.5,0.5);x_1=x0+best_step*search_direction;norm=norm0;k=1;while norm>error && k<80gradient=gradient_my(f,x_1,b);norm=0;for i=1:bnorm=norm+(gradient(i))入2;endnew_direction=-gradient+norm/norm0*search_direction;x=x_1+step_zzh*new_direction;f_step=subs(f,findsym(f),x);best_step=golden_search(f_step,-0.1,0.1)x_2=x_1+best_step*new_directionA=[A;x_2];norm0=norm;search_direction=new_direction;x_1=x_2;k=k+1;end[w。