目标函数极值求解的几种方法题目:()()2221122min -+-x x,取初始点()()Tx 3,11=,分别用最速下降法,牛顿法,共轭梯度法编程实现。
一维搜索法:迭代下降算法大都具有一个共同点,这就是得到点()k x 后需要按某种规则确定一个方向()k d ,再从()k x 出发,沿方向()k d 在直线(或射线)上求目标函数的极小点,从而得到()k x 的后继点()1+k x ,重复以上做法,直至求得问题的解,这里所谓求目标函数在直线上的极小点,称为一维搜索。
一维搜索的方法很多,归纳起来大体可以分为两类,一类是试探法:采用这类方法,需要按某种方式找试探点,通过一系列的试探点来确定极小点。
另一类是函数逼近法或插值法:这类方法是用某种较简单的曲线逼近本来的函数曲线,通过求逼近函数的极小点来估计目标函数的极小点。
本文采用的是第一类试探法中的黄金分割法。
原理书上有详细叙述,在这里介绍一下实现过程:⑴ 置初始区间[11,b a ]及精度要求L>0,计算试探点1λ和1μ,计算函数值()1λf 和()1μf ,计算公式是:()1111382.0a b a -+=λ,()1111618.0a b a -+=μ。
令k=1。
⑵ 若L a b k k <-则停止计算。
否则,当()K f λ>()k f μ时,转步骤⑶;当()K f λ≤()k f μ时,转步骤⑷ 。
⑶ 置k k a λ=+1,k k b b =+1,k k μλ=+1,()1111618.0++++-+=k k k k a b a μ,计算函数值()1+k f μ,转⑸。
⑷ 置k k a a =+1,k k b μ=+1,k k μμ=+1,()1111382.0++++-+=k k k k a b a λ,计算函数值()1+k f λ,转⑸。
⑸ 置k=k+1返回步骤 ⑵。
1. 最速下降法实现原理描述:在求目标函数极小值问题时,总希望从一点出发,选择一个目标函数值下降最快的方向,以利于尽快达到极小点,正是基于这样一种愿望提出的最速下降法,并且经过一系列理论推导研究可知,负梯度方向为最速下降方向。
最速下降法的迭代公式是()()()k k k k d x x λ+=+1,其中()k d 是从()k x 出发的搜索方向,这里取在点()k x 处最速下降方向,即()()k k x f d -∇=。
k λ是从()k x 出发沿方向()k d 进行的一维搜索步长,满足()()()()()()k k k k k d x f d x f λλλ+=+≥0min 。
实现步骤如下:⑴ 给定初点()n R x ∈1 ,允许误差0>ε,置k=1。
⑵ 计算搜索方向()()k k x f d -∇=。
⑶ 若()ε≤k d ,则停止计算;否则,从()k x 出发,沿方向()k d 进行的一维搜索,求k λ,使()()()()()()k k k k k d x f d x f λλλ+=+≥0min 。
⑷ ()()()k k k k d x x λ+=+1,置k=k+1返回步骤 ⑵。
2. 拟牛顿法基本思想是用不包括二阶导数的矩阵近似牛顿法中的Hesse 矩阵的逆矩阵,因构造近似矩阵的方法不同,因而出现了不同的拟牛顿法。
牛顿法迭代公式:()()()k k k k d x x λ+=+1,()k d 是在点()k x 处的牛顿方向,()()()()()k k k x f x f d ∇-∇=-12,k λ是从()k x 出发沿牛顿方向()k d 进行搜索的最优步长。
用不包括二阶导数的矩阵k H 近似取代牛顿法中的Hesse 矩阵的逆矩阵()()12-∇k x f ,1+k H 需满足拟牛顿条件。
实现步骤:⑴ 给定初点()1x ,允许误差0>ε。
⑵ 置n I H =1(单位矩阵),计算出在()1x 处的梯度()()11x f g ∇=,置k=1。
⑶ 令()k k k g H d -=。
⑷ 从()k x 出发沿方向()k d 搜索,求步长k λ,使它满足()()()()()()k k k k k d x f d x f λλλ+=+≥0min ,令()()()k k k k d x x λ+=+1。
⑸ 检验是否满足收敛标准,若()()ε≤+1k y f ,则停止迭代,得到点()1+-=k x x ,否则进行步骤⑹。
⑹ 若k=n ,令()()11+=k x x ,返回⑵;否则进行步骤⑺。
⑺令()()11++∇=k k x f g ,()()()k k k x x p -=+1,()kk k g g q -=+1,()()()()()()()()k k Tk kT k k k k T k T k k k k q H q H q q H q p p p H H -+=+1,置k=k+1 。
返回⑶。
3. 共轭梯度法若()()()k d d d ,,,21 是n R 中k 个方向,它们两两关于A 共轭,即满足()()k j i j i Ad d j T i ,,1,;,0 =≠=,称这组方向为A 的k 个共轭方向。
共轭梯度法的基本思想是把共轭性与最速下降法相结合,利用已知点处的梯度构造一组共轭方向,并沿这组方向进行搜索,求出目标函数的极小点,根据共轭方向的基本性质这种方法具有二次终止性。
实现步骤如下:⑴ 给定初点()1x ,允许误差0>ε,置 ()()11x y =,()()()11y f d -∇=,k=j=1。
⑵ 若()()ε≤j y f ,则停止计算;否则,作一维搜索,求j λ,满足 ()()()()()()j j j j j d y f d y f λλλ+=+≥0min ,令()()()j j j j d y y λ+=+1。
⑶ 若n j <,则进行步骤⑷,否则进行步骤⑸⑷ 令()()()()j j j j d y f d β+-∇=++11,其中()()()()221jj j y f y f ∇∇=+β,置j=j+1,转⑵。
⑸ 令()()11++=n k y x ,()()11+=k x y ,()()()11y f d -∇=,置j=1,k=k+1,转⑵ 。
4. 实验结果1=。
迭用以上三种方法通过Matlab编程得到实验数据。
初始值()()Tx3,1代精度sum(abs(x1-x).^2)<1e-4。
实验结果分析:由上表格可以看到最速下降法需要四次迭代实现所要求的精度,拟牛顿法和共轭梯度法需要三次。
程序:%精确一维搜索法的子函数,0.618(黄金分割)法,gold.m%输入的变量x为初始迭代点是二维的向量,d为初始迭代方向是二维的向量%输出变量是在[0,10]区间上使函数取得极小值点的步长因子function alfa=gold(x,d)a=0;b=10;tao=0.618;lanmda=a+(1-tao)*(b-a);mu=a+tao*(b-a);alfa=lanmda;%初始化f=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;%目标函数m=f;alfa=mu;n=f;while 1if m>nif abs(lanmda-b)<1e-4alfa=mu; returnelsea=lanmda; lanmda=mu; m=n;mu=a+tao*(b-a); alfa=mu;n=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;endelseif abs(mu-a)<1e-4alfa=lanmda; returnelseb=mu; mu=lanmda; n=m;lanmda=a+(1-tao)*(b-a); alfa=lanmda;m=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;endendend%梯度子函数,tidu.m,输入的变量为二维的向量,返回梯度在x处的数值向量function g=tidu(x)%待求解的函数f=(x(1)-2)^2+2*(x(2)-1)^2;%求函数的梯度表达式g=[2*(x(1)-2) 4*(x(2)-1)];x1=x(1); x2=x(2);%最速下降法极小化函数的通用子函数zuisu.m%输入变量为初始的迭代点,输出变量为极小值点function x0=zuisu(x)%判断梯度范数是否满足计算精度1e-4的要求.是,标志变量设为1,输出结果; %否,标志变量设为0if sum(abs(tidu(x)).^2)<1e-4flag=1; x0=x;elseflag=0;end%循环求解函数的极小点while flag==0d=-tidu(x); a=gold(x,d); x=x+a*d%判断梯度范数是否满足计算精度的要求.是,标志变量设为1,输出结果;%否,标志变量设为0,继续迭代if sum(abs(tidu(x)).^2)<1e-4flag=1; x0=x;elseflag=0;endEnd%拟牛顿法极小化函数的通用子函数,gonge.m%输入变量为初始的迭代点,输出变量为极小值点function x0=ninewton(x)%判断梯度范数是否满足计算精度的要求.是,标志变量设为1,输出结果;%否,标志变量设为0,继续迭代if sum(abs(tidu(x)).^2)<1e-4flag=1; x0=x;elseflag=0;end%初始的H矩阵为单位矩阵h0=eye(2);%循环求解函数的极小点while flag==0%计算新的迭代方向d=-h0*tidu(x)'; a=gold(x,d);x1=(x'+a*h0*d)'; s=x1-x;y=tidu(x1)-tidu(x); v=s*y';%校正H矩阵h0=(eye(2)-s'*y./v)*h0*(eye(2)-y'*s./v)+s'*s./v;%判断下一次和上一次迭代点之差是否满足计算精度的要求.是,标志变量设为1,输出结果;否,标志变量设为0,继续迭代if sum(abs(x-x1).^2)<1e-4flag=1; x0=x;elseflag=0;endx=x1end%共轭剃度法极小化函数的通用子函数,gonge.m%输入变量为初始的迭代点,输出变量为极小值点function x0=gonge(x)%判断梯度范数是否满足计算精度的要求.是,标志变量设为1,输出结果;%否,标志变量设为0,继续迭代if sum(abs(tidu(x)).^2)<1e-4flag=1; x0=x;elseflag=0;end%第一次的迭代方法为负梯度方向d1=-tidu(x);a=gold(x,d1);x1=x+a*d1;%循环求解函数的极小点while flag==0g1=tidu(x); g2=tidu(x1);%利用FR公式求解系数batabata=(g2*g2')/(g1*g1');d2=-g2+bata*d1;a=gold(x1,d2);x=x1;x1=x+a*d2%判断下一次和上一次迭代点之差是否满足计算精度的要求.是,标志变量设为1,输出结果;否,标志变量设为0,继续迭代if sum(abs(x1-x).^2)<1e-4flag=1; x0=x1;elseflag=0;endend。