拟牛顿法

拟牛顿法

牛顿法的收敛速度虽然较快，但要求海森矩阵要可逆，要计算二阶导数和逆矩阵，就加大了就算机计算量。为了克服牛顿法的缺点，同时保持较快收敛速度的优点，就产生了拟牛顿法。拟牛顿法是牛顿法的直接推广，通过在试探点附近的二次逼近引进牛顿条件来确定线搜索方向，它主要有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 的正定性，在下面算法步骤中迭代一定次数后，重置初始点和迭代矩阵再