共轭梯度法原理与实现

所用例子：
求解二次目标函数极小点。设

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\qquad

G(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次或者不下降时就再次插入负梯度方向作为搜索方向，从新开始共轭梯度算法。下面的代码就采用了这种思想。

共轭梯度算法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

性能展示

与拟牛顿法http://blog.csdn.net/u012176591/article/details/46225289 对比，发现共轭梯度法还是挺挫的，需要的迭代次数很多，超过一半的样本的迭代次数超过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})

参考文献：

• http://en.wikipedia.org/wiki/Conjugate_gradient_method
• An Introduction to the Conjugate Gradient Method