最优化方法-共轭梯度法
1.简介
共轭梯度法最初由Hesteness和Stiefel于1952年为求解线性方程组而提出的。其基本思想是把共轭性与最速下降方法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜素，求出目标函数的极小点。根据共轭方向基本性质，这种方法具有二次终止性。
对于二次凸函数的共轭梯度法：
m
i
n
f
(
x
)
=
1
2
X
T
A
X
+
B
T
X
+
C
minf(x) = \frac{1}{2}X^TAX+B^TX+C
minf(x)=21XTAX+BTX+C
2.实例
用共轭梯度法求二次函数的极小值与极小点，设初值为[1,1]，迭代精度为0.001
f
(
x
1
，
x
2
)
=
x
1
2
+
2
x
2
2
−
4
x
1
−
2
x
1
x
2
f(x_1，x_2) = x^2_1+2x^2_2-4x_1-2x_1x_2
f(x1，x2)=x12+2x22−4x1−2x1x2
A.计算偏导
▽
f
(
x
0
)
=
[
2
x
1
−
2
x
2
−
4
4
x
2
−
2
x
1
]
=
[
−
4
2
]
\bigtriangledown f(x_0)=\quad \begin{bmatrix} 2x_1-2x_2-4\\4x_2-2x_1\end{bmatrix} \quad=\quad \begin{bmatrix} -4\\2\end{bmatrix} \quad
▽f(x0)=[2x1−2x2−44x2−2x1]=[−42]
B.搜索最佳步长
x
i
=
x
0
+
a
0
▽
f
(
x
0
)
=
=
[
1
1
]
+
a
0
[
−
4
2
]
=
[
1
+
4
a
0
1
−
2
a
0
]
f
(
x
i
)
=
min
a
f
(
x
i
)
=
min
a
(
(
1
+
4
a
)
2
+
2
(
1
−
2
a
)
2
−
4
(
1
+
4
a
)
−
2
(
1
+
4
a
)
(
1
−
2
a
)
)
=
min
a
(
40
a
2
−
20
a
−
3
)
a
0
=
0.25
,
x
1
=
[
2
0.5
]
,
▽
f
(
x
1
)
=
[
−
1
−
2
]
x_i=x_0+a_0\bigtriangledown f(x_0)=\quad=\quad \begin{bmatrix} 1\\1\end{bmatrix} \quad+a_0 \quad \begin{bmatrix} -4\\2\end{bmatrix} \quad= \quad \begin{bmatrix} 1+4a_0\\1-2a_0\end{bmatrix} \quad \\f(x_i)=\min\limits_{a}f(x_i) \\=\min\limits_{a}((1+4a)^2+2(1-2a)^2-4(1+4a)-2(1+4a)(1-2a)) \\=\min\limits_{a}(40a^2-20a-3) \\ a_0=0.25,x_1=\quad \begin{bmatrix} 2\\0.5\end{bmatrix} \quad,\bigtriangledown f(x_1)= \begin{bmatrix} -1\\-2\end{bmatrix} \quad
xi=x0+a0▽f(x0)==[11]+a0[−42]=[1+4a01−2a0]f(xi)=aminf(xi)=amin((1+4a)2+2(1−2a)2−4(1+4a)−2(1+4a)(1−2a))=amin(40a2−20a−3)a0=0.25,x1=[20.5],▽f(x1)=[−1−2]
C.第二次迭代
b
0
=
∥
▽
f
(
x
1
)
∥
2
∥
▽
f
(
x
0
)
∥
2
=
0.25
s
1
=
−
▽
f
(
x
1
)
−
b
0
▽
f
(
x
0
)
=
[
2
1.5
]
x
i
+
1
=
x
i
+
a
▽
f
(
x
i
)
=
=
[
2
0.5
]
+
a
[
2
1.5
]
=
[
2
+
2
a
0.5
+
1.5
a
]
a
=
1.25
,
x
2
=
[
4
2
]
,
▽
f
(
x
2
)
=
[
0
0
]
∥
▽
f
(
x
2
)
∥
2
=
0
<
ξ
收敛
b_0= \frac{\parallel \bigtriangledown f(x_1)\parallel^2}{\parallel \bigtriangledown f(x_0)\parallel^2}=0.25 \\ s1=- \bigtriangledown f(x_1)-b_0\bigtriangledown f(x_0)=\quad \begin{bmatrix} 2\\1.5\end{bmatrix} \quad \\ x_{i+1}=x_i+a\bigtriangledown f(x_i)=\quad=\quad \begin{bmatrix} 2\\0.5\end{bmatrix} \quad+a \quad \begin{bmatrix} 2\\1.5\end{bmatrix} \quad= \quad \begin{bmatrix} 2+2a\\0.5+1.5a\end{bmatrix} \quad \\ a=1.25,x_2=\quad \begin{bmatrix} 4\\2\end{bmatrix} \quad,\bigtriangledown f(x_2)= \begin{bmatrix} 0\\0\end{bmatrix} \quad \\ \parallel \bigtriangledown f(x_2)\parallel^2=0<\xi收敛
b0=∥▽f(x0)∥2∥▽f(x1)∥2=0.25s1=−▽f(x1)−b0▽f(x0)=[21.5]xi+1=xi+a▽f(xi)==[20.5]+a[21.5]=[2+2a0.5+1.5a]a=1.25,x2=[42],▽f(x2)=[00]∥▽f(x2)∥2=0<ξ收敛
3.Matlab代码
二元问题
clc
clear
close all;
syms xi yi;
func = xi^2+2*yi^2-4*xi-2*xi*yi;%创建符号表达式f
f=conjugate_grad_2d(func,[1,1],0.001);
function f = conjugate_grad_2d(func,x0,t)
%用共轭梯度法求已知函数f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2的极值点
%已知初始点坐标：x0
%已知收敛精度：t
%求的已知函数的极值：f
x=x0;
syms xi yi a;%定义自变量，步长为符号变量
% func = xi^2+2*yi^2-4*xi-2*xi*yi;%创建符号表达式f
fx=diff(func,xi);%求表达式的偏导数
fy=diff(func,yi);
dxv=subs(fx,{xi,yi},x);%代入初值计算当前点做多对xi的一阶偏导实值
dyv=subs(fy,{xi,yi},x);
fi=[dxv,dyv];%初始点梯度向量
count=0;%搜索次数初始值
while double(sqrt(dxv^2+dyv^2))>t%搜素精度不满足已知条件
if count<=0 %第一次搜索的方向为负梯度方向
ds=-fi;
else
ds=s1;%第二次搜索方向为-fii+d*s;
end
x=x+a*ds;%进行一次搜索后的点坐标
func_a=subs(func,{xi,yi},x);%构造一元函数φ(a),计算步长
fda=diff(func_a);%对函数求导
dav=solve(fda);%得到最佳步长a
if dav~=0
if length(dav)>1
ai=double(dav(1));
else
ai=double(dav);
end
else
break;%若点a=0,则直接跳出循环，此点为极值
end
%计算下一步梯度
x=subs(x,a,ai);%下一步点坐标
dfxiv=subs(fx,{xi,yi},x);%代入初值计算当前点做多对xi的一阶偏导实值
dfyiv=subs(fy,{xi,yi},x);
dfii=[dfxiv,dfyiv];%下一点梯度向量
b0=(dfxiv^2+dfyiv^2)/(dxv^2+dyv^2);%
s1=-dfii+b0*ds;%下一点搜索的方向向量
count=count+1;%搜索次数+1
dxv=dfxiv;
dyv=dfyiv;
end
x,f=subs(func,{xi,yi},x),count%输出极值点，极小值以及搜索次数
三元问题
clc
clear
close all;
% syms xi yi;
% func = xi^2+2*yi^2-4*xi-2*xi*yi;%创建符号表达式f
% f=conjugate_grad_2d(func,[1,1],0.001);
N=3;
X = sym('x', [1 N])';
% A = rand(N,N);
% B = rand(1,N);
% C = rand(1,N);
% func = X'*A*X +B*X + C;
syms x1 x2 x3;
func = x1^2+2*x2^2+x3^2;%创建符号表达式f-4*x1-2*x1*x2
f=conjugate_grad_2d_new(func,[1,1,1],0.000001);
function f = conjugate_grad_2d_new(func,x0,t)
%用共轭梯度法求已知函数f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2的极值点
%已知初始点坐标：x0
%已知收敛精度：t
%求的已知函数的极值：f
x=x0;
syms x1 x2 x3 a;%定义自变量，步长为符号变量
% func = x1^2+2*x2^2-4*x1-2*x1*x2;%创建符号表达式f
fx=diff(func,x1);%求表达式的偏导数
fx2=diff(func,x2);
fx3=diff(func,x2);
dxv=subs(fx,{x1,x2,x3},x);%代入初值计算当前点做多对x1的一阶偏导实值
dx2=subs(fx2,{x1,x2,x3},x);
dx3=subs(fx3,{x1,x2,x3},x);
fi=[dxv,dx2,dx2];%初始点梯度向量
count=0;%搜索次数初始值
while double(sqrt(dxv^2+dx2^2 + dx3^2))>t%搜素精度不满足已知条件
if count<=0 %第一次搜索的方向为负梯度方向
ds=-fi;
else
ds=s1;%第二次搜索方向为-fii+d*s;
end
%计算当前x的偏导值
dxv=subs(fx,{x1,x2,x3},x);%代入初值计算当前点做多对x1的一阶偏导实值
dx2=subs(fx2,{x1,x2,x3},x);
dx3=subs(fx3,{x1,x2,x3},x);
x=x+a*ds;%进行一次搜索后的点坐标
func_a=subs(func,{x1,x2,x3},x);%构造一元函数φ(a),计算步长
fda=diff(func_a);%对函数求导
dav=solve(fda);%得到最佳步长a
if dav~=0
if length(dav)>1
ai=double(dav(1));
else
ai=double(dav);
end
else
break;%若点a=0,则直接跳出循环，此点为极值
end
%计算下一步梯度
x=subs(x,a,ai);%下一步点坐标
dfx1v=subs(fx,{x1,x2, x3},x);%代入初值计算当前点做多对x1的一阶偏导实值
dfx2v=subs(fx2,{x1,x2, x3},x);
dfx3v=subs(fx3,{x1,x2, x3},x);
dfii=[dfx1v,dfx2v,dfx3v];%下一点梯度向量
b0=(dfx1v^2+dfx2v^2+dfx3v^2)/(dxv^2+dx2^2+dx3^2);%
s1=-dfii+b0*ds;%下一点搜索的方向向量
count=count+1;%搜索次数+1
dxv=dfx1v;
dx2=dfx2v;
dx3=dfx3v;
end
f=subs(func,{x1,x2,x3},x)
double(x),double(f),count%输出极值点，极小值以及搜索次数
%x:1.0e-07 *-0.2878 0.1271 0.1271
%f:1.3134e-15
%c:7
4.非线性共轭梯度
当目标函数是高于二次的连续函数(即目标函数的梯度存在)时，其对应的解方程是非线性方程，非线性问题的目标函数可能存在局部极值，并且破坏了二次截止性，共轭梯度法需要在两个方面加以改进后，仍然可以用于实际的反演计算，但共轭梯度法不能确保收敛到全局极值。
(1)首先是共轭梯度法不能在n维空间内依靠n步搜索到达极值点，需要重启共轭梯度法，继续迭代，以完成搜索极值点的工作。
(2)在目标函数复杂，在计算时，由于需要局部线性化，需计算Hessian矩阵A，且计算工作量比较大，矩阵A也有可能是病态的。
参考文献：
https://www.cnblogs.com/walccott/p/4956966.html
https://wenku.baidu.com/view/613e7216336c1eb91a375da4.html
https://wenku.baidu.com/view/19b537e164ce0508763231126edb6f1afe007173.html
扫一扫
854
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
