作业3-牛顿法

最优化方法第三次作业题目：分别利用最速下降法、阻尼牛顿法、修正牛顿法求解无约束优化问题2122211min2xxxxfRx。该问题有精确解0,1,1**xfxT。初始点分别取50010.0,0,1,1kTTxfxx精确度。比较三种方法的收敛速度。三种方法：1、最速下降法最速下降法是用负梯度方向dk=−∇f(xk)为搜索方向的。步骤：步0选取初始点x0∈Rn,容许误差0≤ε≪1.令k:=1.步1计算gk=∇f(xk).若‖gk‖≤ε,停算,输出xk作为近似最优解.步2取方向dk=−gk.步3由线搜索技术确定步长因子αk.（精确搜索）步4令xk+1:=xk+αkdk,k:=k+1,转步1.程序：function[x,val,k]=grad(fun,gfun,x0)%功能:用最速下降法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun分别是目标函数和梯度%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=5000;%最大迭代次数rho=0.5;sigma=0.4;k=0;epsilon=1e-5;while(kmaxk)g=feval(gfun,x0);%计算梯度d=-g;%计算搜索方向if(norm(d)epsilon),break;endm=0;mk=0;while(m20)%Armijo搜索if(feval(fun,x0+rho^m*d)feval(fun,x0)+sigma*rho^m*g’*d)mk=m;break;endm=m+1;endx0=x0+rho^mk*d;k=k+1;endx=x0;val=feval(fun,x0);解:首先建立两个分别计算目标函数和梯度的m文件:functionf=fun(x)f=(x(1)^2-x(2))^2+(x(1)-1)^2;functiong=gfun(x)g=[4*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-2*(x(1)^2-x(2))]';最速下降法的数值结果.初始点(x0)迭代次数(k)目标函数值(f(xk))(0,0)T1408.5333e-011(−1,−1)T1458.7013e-011由此可以看出，最速下降法的收敛速度是比较慢的。2、阻尼牛顿法步骤：步0给定终止误差值0≤ε≪1,δ∈(0,1),σ∈(0,0.5).初始点x0∈Rn.令k:0.步1计算gk=∇f(xk).若‖gk‖≤ε,停算,输出x*≈xk.步2计算Gk=∇2f(xk),并求解线性方程组得解dk:Gkd=−gk步3记mk是满足下列不等式的最小非负整数m:f(xk+δmdk)≤f(xk)+σδmgTkdk.步4令xk+1=xk+αkdk,k:=k+1,转步1.程序：function[x,val,k]=dampnm(fun,gfun,Hess,x0)%功能:用阻尼牛顿法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun,Hess分别是求%目标函数值,梯度,Hesse阵的函数%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=100;%给出最大迭代次数rho=0.55;sigma=0.4;k=0;epsilon=1e-5;while(kmaxk)gk=feval(gfun,x0);%计算梯度Gk=feval(Hess,x0);%计算Hesse阵dk=-Gk\gk;%解方程组Gk*dk=-gk,计算搜索方向if(norm(gk)epsilon),break;end%检验终止准则m=0;mk=0;while(m20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;endx0=x0+rho^mk*dk;k=k+1;endx=x0;val=feval(fun,x);解:建立的两个计算目标函数和梯度的m文件之外,还需建立求Hesse阵的m文件:functionHe=Hess(x)He=[12*x(1)^2-4*x(2)+2,-4*x(1);-4*x(1),2];阻尼牛顿法的数值结果.初始点(x0)迭代次数(k)目标函数值(f(xk))(0,0)T51.6675e-012(−1,−1)T78.6147e-013由此看出，阻尼牛顿法的收敛速度是比较令人满意的.3、修正牛顿法步骤：步0给定参数δ∈(0,1),τ∈[0,1],σ∈(0,0.5),终止误差0≤ε≪1.初始点x0∈Rn.令k:=0.步1计算gk=∇f(xk),µk=‖gk‖1+τ.若‖gk‖≤ε,停算,输出xk作为近似极小点.步2计算Hesee阵Gk=∇2f(xk).解线性方程组(Gk+µkI)d=−gk,得解dk.步3令mk是满足下列不等式的最小非负整数m:f(xk+δmdk)≤f(xk)+σδmgTkdk.令αk=δmk,xk+1:=xk+αkdk.步4令k:=k+1,转步1.程序：function[x,val,k]=revisenm(fun,gfun,Hess,x0)%功能:用修正牛顿法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun,Hess分别是求%目标函数值,梯度,Hesse阵的函数%输出:x,val分别是近似最优点和最优值,k是迭代次数.n=length(x0);maxk=150;rho=0.55;sigma=0.4;tau=0.0;k=0;epsilon=1e-5;while(k¡maxk)gk=feval(gfun,x0);%计算梯度muk=norm(gk)^(1+tau);Gk=feval(Hess,x0);%计算Hesse阵Ak=Gk+muk*eye(n);dk=-Ak“gk;%解方程组Gk*dk=-gk,计算搜索方向if(norm(gk)epsilon),break;end%检验终止准则m=0;mk=0;while(m20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;endx0=x0+rho^mk*dk;k=k+1;endx=x0;val=feval(fun,x);解:修正牛顿法的数值结果.初始点(x0)迭代次数(k)目标函数值(f(xk))(0,0)T72.7639e-013(−1,−1)T114.1208e-017可以看出,修正牛顿法的收敛速度不如阻尼牛顿法快,因为矩阵Ak=Gk+µkI只是Hesee阵Gk的一个近似.结论：收敛速度最快的是阻尼牛顿法，修正牛顿法其次，最速下降法最慢。