FR共轭梯度法与拟牛顿法计算机实现及仿真.docx
《FR共轭梯度法与拟牛顿法计算机实现及仿真.docx》由会员分享，可在线阅读，更多相关《FR共轭梯度法与拟牛顿法计算机实现及仿真.docx（15页珍藏版）》请在冰点文库上搜索。
FR共轭梯度法与拟牛顿法计算机实现及仿真
FR共轭梯度法与拟牛顿法计算机实现及仿真
信计0701作者：
翟铮学号：
5
摘要：
最速下降法和牛顿法是最基本的无约束最优化方法，但由于最速下降法收敛速度慢而牛顿法虽然收敛速度比较快，但需要计算目标函数的Hesse矩阵及其逆矩阵，计算量比较大。
我们这里采用采用一种无需计算二阶导数，但收敛速度比较快的一种算法，即FR共轭梯度。
关键词:
FR共轭梯度法、拟牛顿法、黄金分割法、最优一维线性搜索、收敛性
一、FR共轭梯度法与拟牛顿法
1、共轭梯度法
在共轭梯度方向法中，如果取初始的搜索方向
，而以以下各共轭方向
由第k次迭代点的负梯度
与已经得到的共轭方向
的线性组合来确定，这样就构造了一种具体的共轭方向法。
因为每一个共轭方向都依赖于迭代点处的负梯度，所以称之为共轭梯度法（conjugategradientmethod）
给定初始点
，取初始搜索方向
，从
出发，沿
进行一维搜索，得到迭代点
，一下按
来构造搜索方向，
的选取应该使得
和
是共轭的。
因为
且要使
和
是共轭的，应有
，所以
于是得到n个搜索方向
如下：
Fletcher_Reeves公式为
2、拟牛顿法：
牛顿法成功的关键是利用了Hesse矩阵提供的曲率信息，但计算Hesse矩阵工作量大，并且有的目标函数的Hesse矩阵很难计算，甚至不好求出。
针对这一问题，拟牛顿法比牛顿法更为有效。
这类算法仅利用目标函数值和一阶导数的信息，构造出目标函数的曲率近似，使方法具有类似牛顿法的收敛速度快的优点。
一般的拟牛顿法算法为：
给出初始点
如果
解
得搜索方向
（或计算
）
由线性搜索得步长因子
校正
使得拟牛顿条件成立。
本次实验采用DFP校正，校正公式为：
二、计算机实现步骤
由于该问题设计的问题比较复杂，所用变量多，为防止变量之间的冲突，在下面的程序中全部用函数来实现。
函数一：
计算函数在任意点处的梯度值。
函数名：
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函数
设定黄金分割法的步长区间为
，搜索精度为0.001，设定共轭梯度法误差为：
0.001
调用函数：
conjugate_gradient（f,[-1.2,1],0.001）;
结论：
从图上可以看到利用共轭梯度法求解过程中，从初始点向收敛点迭代的步长是阶段性的变化的，期间有的阶段步长很小，阶段间的步长比较大。
此方法在搜索在收敛中比较快，有效的避免了最速下降法的锯齿效应。
第三题：
Wood函数
解：
我们用上述的共轭梯度法求解，导致搜索精度低，达不到函数停止迭代的下限，当搜索次数增大时，导致了迭代不收敛，因此我们试图改用拟牛顿法来求解，观察拟牛顿法在处理高维函数时的优异性。
我们用拟牛顿法设定搜索步长区间为[-15,15]，误差精度为：
0.001得到的结果如下：
迭代次数
X
Y
Z
W
1
-3.00000
-1.00000
-3.00000
-1.00000
2
-2.19396
-0.86038
-2.27451
-0.87381
3
-2.44735
-0.57344
-2.05230
-0.54747
4
-3.08965
5.91018
-2.94903
5.31308
5
1.83003
6.50315
0.98570
5.77727
6
-2.58836
7.31350
-2.55296
5.74866
7
-2.07720
4.07977
-2.88577
8.54865
8
-2.08072
3.93519
-2.93265
8.63671
9
-1.42740
1.08110
-3.22100
10.46009
10
-0.80527
0.27410
-3.06216
8.88119
76
1.07314
1.14360
0.91394
0.81642
77
1.06164
1.12478
0.92729
0.84102
78
1.04445
1.09663
0.94733
0.87837
79
1.01963
1.04288
0.96666
0.91740
80
1.01630
1.03885
0.96776
0.92561
81
1.01066
1.03129
0.97072
0.94075
82
0.97768
0.96009
1.01125
1.02000
83
0.96822
0.93931
1.02478
1.04712
84
0.98839
0.97846
1.00975
1.02145
85
1.00025
1.00087
1.00022
1.00027
86
1.00020
1.00042
0.99978
0.99957
87
0.99999
0.99998
1.00000
1.00001
第四题：
Powell奇异函数：
解：
设定黄金分割法的步长区间为：
[-0.1,0.1]，搜索精度为：
0.0001，共轭梯度法的终止点误差范围为
。
迭代终止点为：
[0.0673-0.00670.03320.0334]
调用函数：
conjugate_gradient（f,[3,-1,0,1],0.0001）;可得结果为：
由于函数为四维，我们以下下分别用三个图形来表示收敛效果，分别列出了
图1、X-Y-Z图2、X-Y-W图三Y-Z-W
图1
图2
图3
结论：
从上图可以看到，共轭梯度法搜索过程中，初始步长较大，进行搜索，随着距离收敛点的距离面小，步长逐步变小，向最优解逼近。
共轭梯度法主程序：
functionA=conjugate_gradient（f,x0,error）
[a,b]=size（x0）;
initial_gradient=gradient_my（f,x0,b）;
norm=0;
norm0=0;
symsstep_zzh;
A=[];
fori=1:
b
norm0=norm0+（initial_gradient（i））^2;
end
search_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;
whilenorm>error&&k<80
gradient=gradient_my（f,x_1,b）;
norm=0;
fori=1:
b
norm=norm+（gradient（i））^2;
end
new_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_direction
A=[A;x_2];
norm0=norm;
search_direction=new_direction;
x_1=x_2;
k=k+1;
end
[w,z]=size（A）;
ifz==2
plot（A）;
end
a_return=x;
拟牛顿法主程序：
functionA=quasi_Newton（f,x0,error）
[a,b]=size（x0）;
G0=eye（b）;
initial_gradient=gradient_my（f,x0,b）;
norm0=0;
norm0=initial_gradient*initial_gradient';
symsstep_zzh;
A=[x0];
search_direction=-initial_gradient;
x=x0+step_zzh*search_direction;
f_step=subs（f,findsym（f）,x）;
best_step=golden_search（f_step,-15,15）;
x_1=x0+best_step*search_direction;
A=[A;x_1];
k=1;
whilenorm0>error
ox=x_1-x0;
og=gradient_my（f,x_1,b）-initial_gradient;
G1=G0+（ox'*ox）/（ox*og'）-（G0*og'*og*G0）/（og*G0*og'）;
ifk+1==b
new_direction=-gradient_my（f,x_1,b）;
else
new_direction=-（G1*（gradient_my（f,x_1,b））'）';
end
x=x_1+step_zzh*new_direction;
f_step=subs（f,findsym（f）,x）;
best_step=golden_search（f_step,-15,15）
x_2=x_1+best_step*new_direction
A=[A;x_2];
initial_gradient=gradient_my（f,x_1,b）;
norm0=initial_gradient*initial_gradient';
x0=x_1;x_1=x_2;
G0=G1;
k=k+1;
end
求函数在任意点的梯度向量：
functiongr=gradient_my（f,y,num）
x='';
var='x_';
varible=[];
fori=1:
num
var
（2）=x（i）;
VAR=sym（var）;
varible=[varible,VAR];
end
r=[];
fori=1:
num
r=[r,diff（f,varible（i））];
end
gr=subs（r,findsym（r）,y）;
黄金分割法：
functionbestval=golden_search（f,a,b）
lamda0=a+0.382*（b-a）;
muu0=a+0.618*（b-a）;
f_l=subs（sym（f）,findsym（sym（f））,lamda0）;
f_muu=subs（sym（f）,findsym（sym（f））,muu0）;
iff_l>f_muu
a=lamda0;
else
b=muu0;
end
whileb-a>0.001&&abs（f_l-f_muu）>0.01
lamda=a+0.382*（b-a）;
muu=a+0.618*（b-a）;
f_l=subs（sym（f）,findsym（sym（f））,lamda）;
f_muu=subs（sym（f）,findsym（sym（f））,muu）;
iff_l>f_muu
a=lamda;
else
b=muu;
end
end
bestval=（a+b）/2;
ree=（f_l+f_muu）/2
end
升级会员 