拟牛顿法
牛顿法的收敛速度虽然较快，但要求海森矩阵要可逆，要计算二阶导数和逆矩阵，就加大了就算机计算量。为了克服牛顿法的缺点，同时保持较快收敛速度的优点，就产生了拟牛顿法。拟牛顿法是牛顿法的直接推广，通过在试探点附近的二次逼近引进牛顿条件来确定线搜索方向，它主要有DFP 和BFGS 两种形式，拟牛顿法的一般步骤为：
（1） 给定初始点(0)x ，初始对称正定矩阵0H ，(0)
0()g g x =及精度0ε>； （2） 计算搜索方向()
k k k p H g =-；
（3） 作直线搜索(1)
()()(,)k k k x
F x p +=，计算(1)(1)11(),()k k k k f f x g g x ++++==，
(1)()1,k k k k k k S x x y g g ++=-=-
（4） 判断终止准则是否满足；
（5） 令1k k k H H E +=+置1k k =+，转步骤（2）；
不同的拟牛顿法对应不同的k E ，主要介绍DFP 和BFGS 两种拟牛顿法。
1. DFP 法
（1） 算法原理
DFP 算法中的校正公式为：
1k k
k k
T T k k k k
k k T T k
k k S S H y y H H H S y y H y +=+
-
为了保证k H 的正定性，在下面的算法中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。
（2） 算法步骤
1） 给定初始点(0)x ，初始矩阵0n H I =及精度0ε>； 2） 若()
(0)
f x
ε?≤，停止，极小点为(0)x ；否则转步骤3）；
3） 取()
(0)(0)0p H f x =-?，且令0k =； 4） 用一维搜索法求k t ，使得()
()()()0
()min ()k k k k k k t f X
t p f X tp α≥+=+，令
(1)()()k k k x x tp +=+，转步骤5）；
5） (
)
(1)
k f x
ε+?≤，停止，极小值点为(1)k x +；否则转步骤6）；
6） 若1k n +=，令(0)
()n x x =，转步骤3）；否则转步骤7）；
7） 令
()()()()()
()
()()(
)()()()()()(
)()()()
(1)()(1)()1(1)
()(1)()(1)
()
(1)
()
(1)
()
(1)
()
T
k k k k k k T
k k k k T
k k k k k k
T
k k k k k
x x x x H H x
x f x f x H f x
f x f x f x H f x f x H f x f x +++++++++--=+
-?-??-??-?-
?-??-?，
取()
()(1)1k k k p H f x ++=-?，置1k k =+，转步骤4）。
（3） 算法的MATLAB 程序
调用格式：[,min ]min (,0,var,)x f DFP f x eps = 其中，f ：目标函数 0x ：初始点
var ：自变量向量 eps ：精度
x ：目标函数取最小值时的自变量值 min f ：目标函数的最小值
DFP 的MATLAB 程序代码如下：
function [x,minf]=minDFP(f,x0,var,eps) %目标函数:f; %初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:minf; format long; if nargin==3 eps=1.0e-6; end
x0=transpose(x0); n=length(var); syms l; H=eye(n,n);
gradf=jacobian(f,var); v0=Funval(gradf,var,x0); p=-H*transpose(v0); k=0; while 1
v=Funval(gradf,var,x0); tol=norm(v); if tol<=eps x=x0;
break; end
y=x0+l*p;
yf=Funval(f,var,y); [a,b]=minJT(yf,0,0.1); xm=minHJ(yf,a,b); x1=x0+xm*p;
vk=Funval(gradf,var,x1); tol=norm(vk); if tol<=eps x=x1; break; end
if k+1==n x0=x1; continue; else
dx=x1-x0; dgf=vk-v;
dgf=transpose(dgf); dxT=transpose(dx); dgfT=transpose(dgf); mdx=dx*dxT; mdgf=dgf*dgfT;
fz=H*(dgf*(dgfT*H));
H=H+mdx/(dxT*dgf)-inv(dgfT*(H*dgf))*fz; p=-H*transpose(vk); k=k+1; x0=x1; end end
minf=Funval(f,var,x); format short;
2. BFGS 法
（1） 算法原理
BFGS 算法中的校正公式为
()()()()()()()()()()()1()()()()()()()()()()
11
T T
k k k k k k k T T k k k k T T k k k k k k T
k k S S y H y H H S y S y S y H H y S S y +??
??=++????
??
-
+????
为了保证k H 的正定性，在下面算法步骤中迭代一定次数后，重置初始点和迭代矩阵再
进行迭代。
（2） 算法步骤
1） 给定初始点(0)x ，初始矩阵0n H I =及精度0ε>； 2） 若()
(0)
f x
ε?≤，停止，极小点为(0)x ；否则转步骤3）；
3） 取()
(0)(0)0p H f x =-?，且令0k =； 4） 用一维搜索法求k t ，使得()
()()()0
()min ()k k k k k k t f X
t p f X tp α≥+=+，令
(1)()()k k k x x tp +=+，转步骤5）；
5） (
)
(1)
k f x
ε+?≤，停止，极小值点为(1)k x +；否则转步骤6）；
6） 若1k n +=，令(0)()n x x =，转步骤3）；否则转步骤7）； 7） 令
()()()()()()()()()()()1()()()()()()()()()()
11
T T
k k k k k k k T T k k k k T T k k k k k k T
k k S S y H y H H S y S y S y H H y S S y +??
??=++????
??
-
+????
其中：
()()
()(1)()
()(1)()k k k k k k S x x y f x f x ++=-=?-?，取()
()(1)1k k k p H f x ++=-?，置1k k =+，
转步骤4）。
（3） 算法的MATLAB 程序
调用格式：[,min ]min (,0,var,)x f BFGS f x eps = 其中，f ：目标函数 0x ：初始点
var ：自变量向量 eps ：精度
x ：目标函数取最小值时的自变量值 min f ：目标函数的最小值
BFGS 的MATLAB 程序代码如下：
function [x,minf]=minBFGS(f,x0,var,eps) %目标函数:f; %初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x;
%目标函数的最小值:minf;
format long;
if nargin==3
eps=1.0e-6;
end
x0=transpose(x0);
n=length(var);
syms l;
H=eye(n,n);
gradf=jacobian(f,var);
v0=Funval(gradf,var,x0);
p=-H*transpose(v0);
k=0;
while 1
v=Funval(gradf,var,x0);
tol=norm(v);
if tol<=eps
x=x0;
break;
end
y=x0+l*p;
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);
x1=x0+xm*p;
vk=Funval(gradf,var,x1);
tol=norm(vk);
if tol<=eps
x=x1;
break;
end
if k+1==n
x0=x1;
continue;
else
dx=x1-x0;
dgf=vk-v;
dgf=transpose(dgf);
dxT=transpose(dx);
dgfT=transpose(dgf);
mdx=dx*dxT;
mdgf=dgf*dgfT;
H=H+(1+dgfT*(H*dgf)/(dxT*dgf))*mdx/(dxT*dgf)-(dx*dgfT*H+H*dgf*dxT)/(dxT*dgf);
p=-H*transpose(vk); k=k+1;
x0=x1;
end
end
minf=Funval(f,var,x); format short;