- 共轭方向
- 定义
- 共轭方向的性质
- 共轭方向法
- 算法描述
- 算法的收敛性
- 搜索步长kalpha_k的确定
- 共轭梯度法
- 共轭梯度法的原理
- 共轭梯度算法描述
- 共轭梯度算法Python实现
所用例子：
求解二次目标函数极小点。设
minf(x)=12xTGx+bTx+c
其中
G是n阶对称正定矩阵，
b为一维常向量，c为常数。
1.共轭方向
定义：
设G是n阶对称正定矩阵，若n维向量组d1,d2,⋯,dm(m≤n)满足：
dTiGdj=0,i≠j
则称
d1,d2,⋯,dm为关于
G共轭的。
当G=I时，则上式变为
dTidj=0,i≠j
即向量相互正交。由此可见共轭概念是正交概念的推广，正交概念是共轭概念的特例。
共轭方向的性质
- 若非零向量d1,d2,⋯,dm对于对称正定矩阵G共轭，则这m个向量线性无关。
- 在n维空间中互相共轭的非零向量不超过n个。
- 从任意初始点出发，依次沿n个G的共轭方向d1,d2,⋯,dm进行一维寻优，最多经过n次寻优就可以找到二次函数的极小值点。
2.共轭方向法
算法描述
step 1 : 给定迭代精度0≤ϵ≪1和初始点x0. 计算g0=∇f(x0). 选取初始方向d0，使得dT0g0<0. 令k←0.
step 2 : 若||gk||≤ϵ，停止迭代，输出x∗≈xk
step 3 : 确定搜索步长αk
step 4 : 令xk+1←xk+αkdk，并计算gk+1=∇f(xk+1)
step 5 : 选取dk+1满足如下下降性和共轭性条件：dTk+1gk+1<0,dTk+1Gdi=0,i=0,1,⋯,k
step 6 : k←k+1，转step 2
算法的收敛性
设目标函数为之前定义的f(x)，{xk}是有算法产生的迭代序列，则每一步迭代xk+1都是f(x)在x0和方向d0,d1,⋯,dk所形成的线性流形
Sk={X|X0+∑i=0kαidi,∀αi}
中的极小点。特别地，xn=x∗=−G−1b是目标函数的唯一极小值点。
证明：有xk+1=xkdk=⋯=x0+∑i=0kαidi∈Sk
设任意x∈Sk，存在γi∈R(i=0,1,⋯,k),使得
x=x0+∑i=0kγidi
记x与xk+1的差为 hk+1，有
hk+1=x−xk+1=∑i=0k(γi−αi)di
利用泰勒展开公式，有
f(x)=fk+1+gTk+1hk+1+12hTk+1Ghk+1≥fk+1+gTk+1hk+1=fk+1+∑i=0k(γi−αi)gTk+1di
下面只需证明
gTk+1di=0,∀i=0,1,⋯,k
即可。实际上，因
gj+1−gj=G(xj+1−xj)=αiGdj
故当i≤k时有
gTk+1di=gTi+1di+∑j=i+1k(gj+1−gj)Tdi=gTi+1di+∑j=i+1kαjdTjdi=0
故每一步迭代xk+1都是f(x)在x0和方向d0,d1,⋯,dk所形成的线性流形
Sk={X|X0+∑i=0kαidi,∀αi}
中的极小点。
搜索步长αk的确定
设x∗是目标函数的极小值点，x0为不同于x∗的任意一点，则它们的差向量可以表示为
x∗−x0=∑i=0n−1αidi
将其改写成如下形式
x∗=x0+∑i=0n−1αidi
然后从逐步寻优的角度分析该式，可以把x0看成初值，按照上式进行n次累加后得到的结果。这是一种经过特殊迭代关系的寻优过程，其中经过k次寻优得到的点xk的计算通式可以表示为
xk=x0+∑i=0k−1αidi
可以将视dk和αk 为k+1次迭代的搜索方向和步长。
对向量x∗−x0左乘dTkG，得到
dTkG(x∗−x0)=∑i=0n−1αidTkGdi=αkdTkGdk
进而得到步长αk的表达式
αk=dTkG(x∗−x0)dTkGdk
对向量xk−x0左乘dTkG，得到
dTkG(xk−x0)=∑i=0k−1αidTkGdi=0
从而得到
dTkGxk=dTkGx0
将该等式带入到αk的表达式，得
αk=dTkG(x∗−xk)dTkGdk
对二次目标函数，其在x处的梯度向量为g(x)=Gx+b，所以Gx=g(x)−b，有
G(x^*-x_k)=\left[g(x^*)-b\right]-\left[g(x_k)-b\right]\\ \qquad\quad=g(x^*)-g(x_k)=0-g(x_k)\\ =-g(x_k)\qquad\qquadG(x∗−xk)=[g(x∗)−b]−[g(xk)−b]=g(x∗)−g(xk)=0−g(xk)=−g(xk)
最后得到
αk=−dTkg(xk)dTkGdk
用共轭方向法的思想可以解决前面给出的二次目标函数f(x)=12xTGx+bTx+c的极小值，这等同于求线性方程组Gx=b的解。
3.共轭梯度法
共轭梯度法的原理
在寻优过程中利用当前点xk处的梯度向量gk和前一迭代点xk处的搜索方向dk−1对搜索方向进行如下修正：
dk=−gk+βk−1dk−1
其修正系数βk−1的取值有一个约束条件，即要确保dk与dk−1,dk−2,⋯,d0之间满足关于G的共轭关系。这就是共轭梯度法的基本思想。
修正系数βk−1的取值方法有多个，下面的例子采用的取值公式为
βk−1=gTkgkgTk−1gk−1
.
可以看出共轭梯度法的搜索方向dk的计算只需要梯度向量，不需要矩阵G，可以推广到非二次目标函数的极小值求解，但是这种推广也带来了构造的搜索向量序列{dk}不共轭的问题，后面有提到解决办法。
共轭梯度算法描述
step 1 : 给定迭代精度0≤ϵ≪1和初始点x0. 计算g0=∇f(x0). . 令k←0.
step 2 : 若||gk||≤ϵ，停止迭代，输出x∗≈xk
step 3 : 计算搜索方向 dk
dk={−gk−gk+βk−1dk−1k=0k≥1
step 4 : 利用线搜索方法确定搜索步长 αk
step 5 : 令 令 xk+1←xk+αkdk，并计算 gk+1=∇f(xk+1)
step 6 : k←k+1，转 step 2说明：
通常来说，共轭梯度法的收敛速度比最速下降法快，而且不用像牛顿法那样计算海森矩阵及其逆矩阵。但是随着迭代次数的增加，新构造的共轭方向由于误差（如果目标函数不是二次函数则会造成这种误差）积累会逐渐不精确甚至不下降，可能出现收敛速度极慢的现象。为了避免这种现象，一种有效的改进办法是：
每迭代n<script type="math/tex" id="MathJax-Element-121">n</script>次或者不下降时就再次插入负梯度方向作为搜索方向，从新开始共轭梯度算法。下面的代码就采用了这种思想。共轭梯度算法Python实现
def frcg(fun,gfun,x0): #用FR共轭梯度法求解无约束问题 #x0是初始点，fun和gfun分别是目标函数和梯度 #x,val分别是近似最优点和最优值，k是迭代次数
maxk = 5000
rho = 0.6
sigma = 0.4
k = 0
epsilon = 1e-5
n = np.shape(x0)[0]
itern = 0
while k < maxk:
gk = gfun(x0)
itern += 1
itern %= n
if itern == 1:
dk = -gk
else:
beta = 1.0*np.dot(gk,gk)/np.dot(g0,g0)
dk = -gk + beta*d0
gd = np.dot(gk,dk)
if gd >= 0.0:
dk = -gk
if np.linalg.norm(gk) < epsilon:
break
m = 0
mk = 0
while m < 20:
if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
x0 += rho**mk*dk
g0 = gk
d0 = dk
k += 1
return x0,fun(x0),k性能展示
与拟牛顿法对比，发现共轭梯度法还是挺挫的，需要的迭代次数很多，超过一半的样本的迭代次数超过500（上图没有显示）。
作图代码：
n = 50x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]
iters = []
for i in xrange(np.shape(data)[0]):
rt = frcg(fun,gfun,data[i])
if rt[2] <=200:
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=50)
plt.title(u'共轭梯度法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})
