共轭梯度法
对于任意形式的目标函数f(X),在极值点X附近展开成泰勒级数,且取前三项,有
T1***T2**
f(X)fXfX.XXXX.fX.XX2
*
2**
因在极值点X处fX*0,而fXH(X)为f(X)在X的二阶偏导数矩阵,
*
*
*
即Hessian矩阵,故
1*T**
f(X)fX*XX.H(X).XX 2
对于二次函数来说,若令
2fX*x1
则
2
a,
2fX*x1x2
b,
2fX*x2
2
c
ab
H(X*),而fX*d1—常数
bc
则,得到
1
f(X)d1+x1-x1*
2
2
1
=d1+x1-x1*
2
x1-x*ab1x2-x*
2*
bcx2-x2
ax1-x*bx2-x*
12
x2-x*2bx1-x*cx2-x*
12
1
d1ax1-x1*
2
2bx1-x1*
2
*x2-x*cx-x222
x1*x1*
由上式可知,当XX时,得到目标函数的极小值
*xx22
f(X)fX*d1,当f(X)d2,d2,...时,则有等值线族。
令f(X)d2,代入上式,则有
f(X)d2d1
*
1
ax1-x1*2
2
2bx1-x1*
2
* x2-x*cx-x222
所以目标函数f(X)在X点附近的等值线方程为
*
ax1-x1
2
*
2bx1-x1
*
x2-x*cx-x222
2
d0
式中,d2(d1d2)常数。由极小值点存在的充分必要条件为二阶导数行列式小于零,
即
b2ac0
在解析几何中已经证明,当b2ac0时为椭圆,所以二元函数的等值线在极值点附近是近似的同心椭圆族。
1. 算法原理
共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法,每次搜索方向都是在目标函数梯度的共轭方向,搜索步长通过一维极值算法确定。
T
设二次函数为f(X)CbX
1T
XAX,其中C为常数,b,X为n维列向量,A2
为对称正定矩阵,用共轭梯度法求f(X)的极小点:
共轭梯度法探索的第一步是沿负梯度方向。即X
(k)
点按S
(k)
fX(k)方向找到
X(k1),然后沿着与上一次探索方向S(k)相共轭的方向S(k1)进行探索直达到最小点X*。
令S
(k1)
fX(k1)kS(k)。
(k)
S(k)的一部分即kS(k),加上新的负梯上式的意义就是以原来的负梯度fX
(k1)
度fX,构造S(k1)。
在上式中的选择,应使n维欧氏空间E中的两个非零向量S(k)与S(k1)关于矩阵A共轭。即
(k1)(k)
SAS0
T
kn
(k0,1,2,...n1)
因
f(X)CbTX
若令
1T
XAX,故有f(X)bAX 2
g(k)fX(k)bAX(k) g(k1)fX(k1)bAX(k1)
二式相减,得
g(k1)g(k)A(X(k1)X(k))
设在第k1次迭代中
X(k1)X(k)(k)S(k)
代入上式,得
g(k1)g(k)A(k)S(k),(k0,1,2,...n1)
式中(k)为第k1次迭代的最优步长。
由式S
(k1)T
(k)
AS0
(k0,1,2,...n1)和式
g(k1)g(k)A(k)S(k),(k0,1,2,...n1),得
(k1)(k1)
g(k))0 S(g
T
将式S
(k1)
fX(k1)kS(k)和式g(k1)fX(k1)bAX(k1)代入上式,得
[g(k1)kS(k)]T(g(k1)g(k))0
因为g(k1),g(k),...,g(0)是一正交系,故有[g(k1)]Tg(k)0或[g(k1)]TS(k)0
故上式可简化为
[g(k1)]Tg(k1)(k)[g(k)]Tg(k)0
得
(k)
g[g]g
[g(k)]Tg(k)g(k)
(k1)T
(k1)
(k1)2
2
fX(k1)fX
(k)
22
(k)用一维探索最优化方法确定,即求
minf(X(k)S(k))f(X(k)(k)S(k))
a
或用解析法,使
df(X(k)(k)S(k))
0
d
求得
(k)
。
(k1)
或由式g
g(k)A(k)S(k),(k0,1,2,...n1),得
g(k1)g(k)A(k)S(k)
即
fX(k1)fX(k)(k)AS(k)
又由于进行一维最优化搜索,故有
fX(k1)S(k)0
T
由上两式可求得
(k)
fX(k)S(k)
(k)T(k)SAS
T
如此,即可得到共轭梯度法的一组计算公式为
X(k1)X(k)(k)S(k)
(k)
fX(k)S(k)fX(k)S(k)
(k)T(k)(k)T(k)(k)
SASSHXS
S(k1)fX(k1)kS(k)
TT
(k)
g[g]g
(k)T(k)
[g]gg(k)
(k1)T
(k1)
(k1)2
2
fX
(k1)(k)
fX
22
2. 算法步骤
用共轭梯度法求无约束多维极值问题minf(x),xRn的算法步骤如下: (1) 给定初始点x(0),及精度0;
(0)(0)
(2) 若f(x),停止,极小值点为x,否则转步骤(3);
(3) 取p
(0)
f(x(0)),且置k0;
(k)
(k)
t0
(4) 用一维搜索法求tk,使得f(xktp)minf
k())
xtkp令,(,
x(k1)x(k)tkp(k),转步骤5;
(k1)
),停止,极小值点为x(k1),否则转步骤(6)(5) 若f(x;
(6) 若k1n,令x
(0)
x(n),转步骤(3),否则转步骤(7);
k(7) 令p(k1)f(x(k1))kp(k),
f(x(k1))f(x)
(k)
22
,置kk1,转步骤(4)。
3. 算法的MATLAB实现
在MATLAB中编程实现的共轭梯度法函数为:minGETD
功能:用共轭梯度法求解多维函数的极值。
调用格式:[x,minf]minGETD(f,x0,var,eps) 其中,f:目标函数;
x0:初始点; var:自变量向量;
x:目标函数取最小值时的自变量值; minf:目标函数的最小值。 共轭梯度法的MATLAB程序代码如下: function [x,minf]=minGETD(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;
gradf=jacobian(f,var); %梯度方向 v0=Funval(gradf,var,x0); p=-transpose(v0); k=0; while 1
v=Funval(gradf,var,x0); tol=norm(v); if tol
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
end
if k+1==n x0=x1; continue; else
lamda=dot(vk,vk)/dot(v,v);
p=-transpose(vk)+lamda*p; %共轭方向 k=k+1; x0=x1; end end
minf=Funval(f,var,x); format short; 例:
用共轭梯度法求函数f(t,s)(t3)2s2的极小值,其中初始值取x0(2,6)。 解:在MATLAB命令窗口中输入:
syms t s;
f=(t-3)^2+s^2;
[x,mf]=minGETD(f,[-2 6],[t s])
所得结果为: x =
3.0000 0.0000 mf =
2.0116e-037
22
例:试用共轭梯度法求解目标函数f(X)x1x2x1x210x14x260的极小值。设初
始点X
(0)
[x1(0)
(0)T
x2]00。
T
解:原函数的梯度为
f(X)[
f(X)
x1f(X)TT
]2x1x2102x2x14 x2f(X(0))104
T
第一次探索的方向:
S(0)g(0)104
T
2f(X)x2
1
2f(X)H(X)2
f(X)x2x1
T
2f(X)
x1x221
2f(X)12
2
x2
T
由式
(k)
fX(k)S(k)fX(k)S(k)
,得 (k)T(k)(k)T(k)(k)
SASSHXS
(0)
fX(0)S(0)fX(0)S(0)
(0)T(0)(0)T
HX(0)S(0)SASS
10
104
11640.763157894
2110152
104
124
TT
由此得
X
求X
(1)
X
(0)
S
(0)(0)
010T0.7631578947.631578943.052631576 04
(1)
点处的梯度为
g(1)2.210526305.52631579
T
(k)
由式
g[g]g
[g(k)]Tg(k)g(k)
(k1)T
(k1)
(k1)2
2
fX(k1)fX
(k)
22
求
(0)
:
(0)
第二次探索的方向:
g
(1)2
2
g(0)
35.42659278
0.305401661
116
S(1)g(1)(0)S(0)0.84349036.74792243
T
由式
(k)
fXS
(k)T(k)SAS
(k)
T
(k)
fXS(k)
T求探索步长为:
(k)(k)(k)
SHXS
(k)
T
0.8434903
2.210526305.52631579
6.74792243(1)
a
210.8434903
0.84349036.747922436.74792243
12
35.426592830.[1**********].10825193
由此得
X
(2)
X
(2)
aS
(1)(1)
7.99999995.9999999
T
T
x1(2)(2) x2
g(2)0.00000020.0000002
g(2)
则,极小值点为
2
x1(2)T
(2)7.99999995.9999999 x2