优化方法之正定二次函数的共轭梯度法及其实现（基于Python）

共轭梯度法也是共轭方向法中的一种，但是它减少了梯度方向的搜索量，它直接采取经过一维搜索最小点处的梯度方向作为我们的搜索方向，因而在计算速度上有了一定的提升。如果你对这些优化算法感到困惑，现在你需要明白共轭方向法是基于最速下降法的改进，因为最速下降法在接近最优值时的锯齿现象降低了迭代搜索的速度，共轭法则提升了最速下降法的速度。本节所讲的共轭梯度法则是共轭方向法的进一步改进，一个直观原因是它减少了方向的搜索。同时此算法需要同时计算当前点和下一点并据此得到一些新的中间变元。在你明白这个逻辑 后，请参考以下算法原理和步骤：

结合上述算法步骤，我选取了一个二维函数编写了共轭梯度法的程序进行了验证，相关程序如下：

from matplotlib import pyplot as plt
import numpy as np
import math

#正定二次函数的共轭梯度法，这里共轭梯度法也是共轭方向法的一种，它的搜索方向是利用一维搜索所得极小值点处的梯度生成的
#这个程序指定一个Function二维向量函数以测试该算法，并给出相对应的优化结果。

#定义目标函数和偏导数
def Function(x1,x2):
    value=3*math.pow(x1,2)+0.5*math.pow(x2,2)-x1*x2-2*x1
    return value
def Jac_x1(x1,x2):
    value=6*x1-2-x2
    return value
def Jac_x2(x1,x2):
    value=x2-x1
    return value

#数据初始化
X=[1,2] #初始点k=0
x=[X[0],X[1]] #保留初始点的副本
Hessian=[[6,-1],[-1,1]] #目标函数的二阶偏导矩阵
Jac=[Jac_x1(X[0],X[1]),Jac_x2(X[0],X[1])] #雅可比一阶偏导向量
Direction=[-Jac[0],-Jac[1]]  #初始化搜索方向
#初始化步长
Lambda=-(Jac[0]*Direction[0]+Jac[1]*Direction[1])/((Direction[0]*Hessian[0][0]+Direction[1]*Hessian[1][0])*Direction[0]+
        (Direction[0]*Hessian[0][1]+Direction[1]*Hessian[1][1])*Direction[1]) #计算最优步长lambda

X1=[X[0]+Lambda*Direction[0],X[1]+Lambda*Direction[1]] #计算下一个点
Jac1=[Jac_x1(X1[0],X1[1]),Jac_x2(X1[0],X1[1])]
Beta=(math.pow(Jac1[0],2)+math.pow(Jac1[1],2))/(math.pow(Jac[0],2)+math.pow(Jac[1],2)) #初始化Beta的值

#设置精确度
Epsilon=0.01

result=[]
Minimum=Function(X[0],X[1])
error=[]
count=0

while (math.sqrt(math.pow(Jac[0],2)+math.pow(Jac[1],2))>Epsilon):
    result.append(Function(X[0],X[1]))
    count+=1
    if(Function(X[0],X[1])<Minimum):
        Minimum=Function(X[0],X[1])
        
    #更新各个值
    X[0]=X1[0]
    X[1]=X1[1]
    #更新当前雅可比向量
    Jac=[Jac_x1(X[0],X[1]),Jac_x2(X[0],X[1])] 
    #更新搜索方向
    Direction=[-Jac[0],-Jac[1]] 
    #更新Lambda步长值
    Lambda=-(Jac[0]*Direction[0]+Jac[1]*Direction[1])/((Direction[0]*Hessian[0][0]+Direction[1]*Hessian[1][0])*Direction[0]+
        (Direction[0]*Hessian[0][1]+Direction[1]*Hessian[1][1])*Direction[1])
    #更新当前点的下一点
    X1=[X[0]+Lambda*Direction[0],X[1]+Lambda*Direction[1]]
    Jac1=[Jac_x1(X1[0],X1[1]),Jac_x2(X1[0],X1[1])]
    #更新beta值
    Beta=(math.pow(Jac1[0],2)+math.pow(Jac1[1],2))/(math.pow(Jac[0],2)+math.pow(Jac[1],2))

#存储误差
for i in range(0,count,1):
    error.append(result[i]-Minimum)
#结果分析
print("迭代次数为：{}".format(count))
print("初始点为：（{},{}）时的最优值为：{}".format(x[0],x[1],Minimum))
print("最优点为：（{},{}）".format(X[0],X[1]))
plt.plot(result)
plt.title("The function value in the process of iteration")
plt.xlabel("The number of iteration")
plt.ylabel("The value of the function")
plt.figure()
plt.plot(error)
plt.title("The error in the process of iteration")
plt.xlabel("The number of iteration")
plt.ylabel("The value of error")
    

程序的运行结果见下：

可以看到当我们选择初始迭代点为（1，2）时该算法最终迭代得到的优化值约为-0.4，迭代次数为19次。迭代过程的函数值变化如下图所示：

从上图可见，共轭梯度算法在某一点都可以取得找到对应的下降方向并且最终收敛到一个近似最优值。（读者如果乐意， 你可以直接拷贝本文中的程序并修改程序的精度，初始点甚至函数表达式等，以观察该算法对迭代过程的影响）
该算法的误差变化如下图所示：

可以看到我们的误差也是在逐渐减小的，因此我们的算法效果还是很不错滴，快快行动起来吧！（编程序）