目标函数的几种极值求解方法.doc

目标函数极值求解的几种方法
题目：，取初始点，分别用最速下降法，你牛顿法，共轭梯度法编程实现。
一维搜索法：
迭代下降算法大都具有一个共同点，这就是得到点后需要按某种规则确定一个方向，再从出发，沿方向在直线（或射线）上求目标函数的极小点，从而得到的后继点，重复以上做法，直至求得问题的解，这里所谓求目标函数在直线上的极小点，称为一维搜索。
一维搜索的方法很多，归纳起来大体可以分为两类，一类是试探法：采用这类方法，需要按某种方式找试探点，通过一系列的试探点来确定极小点。另一类是函数逼近法或插值法：这类方法是用某种较简单的曲线逼近本来的函数曲线，通过求逼近函数的极小点来估计目标函数的极小点。本文采用的是第一类试探法中的黄金分割法。原理书上有详细叙述，在这里介绍一下实现过程：
⑴ 置初始区间[]及精度要求L>0,计算试探点和，计算函数值和，计算公式是：，。令k=1。
⑵ 若则停止计算。否则，当>时,转步骤⑶;当时,转步骤⑷ 。
⑶ 置,,,,计算函数值,转⑸。
⑷ 置,,,,计算函数值,转⑸。
⑸ 置k=k+1返回步骤 ⑵。
1. 最速下降法
实现原理描述：在求目标函数极小值问题时,总希望从一点出发,选择一个目标函数值下降最快的方向,以利于尽快达到极小点,正是基于这样一种愿望提出的最速下降法,并且经过一系列理论推导研究可知,负梯度方向为最速下降方向。
最速下降法的迭代公式是，其中是从出发的搜索方向，这里取在点处最速下降方向，即。是从出发沿方向进行的一维搜索步长，满足。
实现步骤如下：
⑴ 给定初点 ，允许误差，置k=1。
⑵ 计算搜索方向。
⑶ 若，则停止计算；否则，从出发，沿方向进行的一维搜索，求，使。
⑷ ，置k=k+1返回步骤 ⑵。
2. 拟牛顿法
基本思想是用不包括二阶导数的矩阵近似牛顿法中的Hesse矩阵的逆矩阵，因构造近似矩阵的方法不同，因而出现了不同的拟牛顿法。
牛顿法迭代公式：，是在点处的牛顿方向，，是从出发沿牛顿方向进行搜索的最优步长。用不包括二阶导数的矩阵近似取代牛顿法中的Hesse矩阵的逆矩阵，需满足拟牛顿条件。
实现步骤：
⑴ 给定初点 ，允许误差。
⑵ 置（单位矩阵），计算出在处的梯度，置k=1。
⑶ 令。
⑷ 从出发沿方向搜索，求步长，使它满足，令。
⑸ 检验是否满足收敛标准，若，则停止迭代，得到点，否则进行步骤⑹。
⑹ 若k=n，令，返回⑵；否则进行步骤⑺。
⑺令，，，，置k=k+1 。返回⑶。

3. 共轭梯度法
若是中k个方向，它们两两关于A共轭，即满足 ，称这组方向为A的k个共轭方向。共轭梯度法的基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出目标函数的极小点，根据共轭方向的基本性质这种方法具有二次终止性。
实现步骤如下：
⑴ 给定初点 ，允许误差，置
，，k=j=1。
⑵ 若，则停止计算；否则，作一维搜索，求，满足
，令。
⑶ 若，则进行步骤⑷，否则进行步骤⑸
⑷ 令，其中，置j=j+1，转⑵。
⑸ 令，，，置j=1，k=k+1，转⑵ 。

4. 实验结果
用以上三种方法通过Matlab编程得到实验数据。初始值 。迭代精度sum(abs(x1-x).^2)<1e-4。

最速下降法
拟牛顿法
共轭梯度法
第一次迭代
结果

1.5151631286
1.5151631286
1.5151631286

0.9393474854
0.939347485
0.9393474854
第二次迭代
结果

1.9730082275
2.0108108072
2.0000076259

1.0538992374
0.9861577108
1.0000419788
第三次迭代
结果

1.9869133934
2.005410162
2.0000038167

0.9983654378
0.9896269240
0.9999998271
第四次迭代
结果

1.9992739761

1.0014531964

实验结果分析：
由上表格可以看到最速下降法需要四次迭代实现所要求的精度，拟牛顿法和共轭梯度法需要三次。

程序：
%精确一维搜索法的子函数,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 1
if m>n
if abs(lanmda-b)<1e-4
alfa=mu; return
else
a=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;
end
else
if abs(mu-a)<1e-4
alfa=lanmda; return
else
b=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;
end
end
end
%梯度子函数，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，输出结果;
%否,标志变量设为0
if sum(abs(tidu(x)).^2)<1e-4
flag=1; x0=x;
else
flag=0;
end
%循环求解函数的极小点
while flag==0
d=-tidu(x); a=gold(x,d); x=x+a*d
%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;
%否,标志变量设为0，继续迭代
if sum(abs(tidu(x)).^2)<1e-4
flag=1; x0=x;
else
flag=0;
end
End

%拟牛顿法极小化函数的通用子函数，gonge.m
%输入变量为初始的迭代点，输出变量为极小值点
function x0=ninewton(x)
%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;
%否,标志变量设为0，继续迭代
if sum(abs(tidu(x)).^2)<1e-4
flag=1； x0=x;
else
flag=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-4
flag=1; x0=x;
else
flag=0;
end
x=x1
end
%共轭剃度法极小化函数的通用子函数，gonge.m
%输入变量为初始的迭代点，输出变量为极小值点
function x0=gonge(x)
%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;
%否,标志变量设为0，继续迭代
if sum(abs(tidu(x)).^2)<1e-4
flag=1; x0=x;
else
flag=0;
end
%第一次的迭代方法为负梯度方向
d1=-tidu(x);a=gold(x,d1);x1=x+a*d1;
%循环求解函数的极小点
while flag==0
g1=tidu(x); g2=tidu(x1);
%利用FR公式求解系数bata
bata=(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-4
flag=1; x0=x1;
else
flag=0;
end
end