共轭梯度法也是解决无约束优化问题的常用迭代算法,它结合了最速下降法矩阵共轭梯度的性质,可以加快算法的迭代过程。且如果初始点选取后的最终优化中不满足精度条件,还可保存上一步得到的迭代点进行再次迭代直到获得较好的优化值。以上过程一般都可以获得较好的迭代点和优化值。该算法简介如下:
根据以上算法过程,我们可以选取目标函数进行测试,以下是测试代码:
(注:读者可以自由地在初始数据修改初始点、精度等参数,以观察和比较不同初始点的迭代过程,迭代最优值的差异等,如果读者希望修改目标函数,则需要在相应的函数定义处进行修改,并注意修改相对应的初始值表达式(如果你改变自变量的个数的话,测试的函数使用的是两个,因此如果你考虑的新函数仍然是二元函数,你则无需修改对应的函数表达式)
import math
import numpy
from matplotlib import pyplot as plt
def Function(x1,x2):
value=1.5*math.pow(x1,2)-x1*x2-2*x1+0.5*math.pow(x2,2)
return value
def Jac_x1(x1,x2):
value=3*x1-x2-2
return value
def Jac_x2(x1,x2):
value=x2-x1
return value
#设置初始迭代值
X=[-3,2] #初始迭代点
x=[X[0],X[1]] #记录初始点
Jac=[Jac_x1(X[0],X[1]),Jac_x2(X[0],X[1])]
Hessian=[[3,-1],[-1,1]]
Direction=[-Jac[0],-Jac[1]]
Lambda=(math.pow(Jac[0],2)+math.pow(Jac[1],2))/(((Jac[0]*Hessian[0][0]+Jac[1]*Hessian[1][0])*Jac[0])+(Jac[0]*Hessian[0][1]+Jac[1]*Hessian[1][1])*Jac[1]) #计算最优迭代步长,利用共轭梯度加最速下降法的思想
#声明存储计算结果的矩阵
result=[]
error=[]
count=0 #迭代次数
Minimum=Function(X[0],X[1]) #记录最小值并初始化为初始点处的值
Epsilon=0.01
while ((math.pow(Jac[0],2)+math.pow(Jac[1],2))>Epsilon):
#记录函数值
result.append(Function(X[0],X[1]))
count+=1
#print(result)
#迭代更新点和相关系数
X[0]=X[0]+Lambda*Direction[0]
X[1]=X[1]+Lambda*Direction[1] #更新点
#更新方向
Direction=[-Jac_x1(X[0],X[1]),-Jac_x2(X[0],X[1])]
#更新Jac矩阵的值
Jac=[-Direction[0],-Direction[1]]
#更新Lambda步长的值
Lambda=(math.pow(Jac[0],2)+math.pow(Jac[1],2))/(((Jac[0]*Hessian[0][0]+Jac[1]*Hessian[1][0])*Jac[0])+(Jac[0]*Hessian[0][1]+Jac[1] *Hessian[1][1])*Jac[1]) #计算最优迭代步长,利用共轭梯度加最速下降法的思想
for i in range (0,count,1):
error.append(result[i]-Minimum) #记录误差函数值
print("初始点选择{}{}".format(x[0],x[1]))
print("最小点是:{},{}".format(X[0],X[1]))
print("最小值是:{}".format(Minimum))
print("该初始点下的迭代值为:{}".format(count))
plt.plot(result)
plt.title("The changes of function value when choosing the initial point:{}{}".format(x[0],x[1]))
plt.xlabel("The number of iteration")
plt.ylabel("The value of function")
plt.show()
plt.figure()
plt.plot(error)
plt.title("The error in the process when choosing initial value:{}{}".format(x[0],x[1]))
plt.xlabel("The number of iteration")
plt.ylabel("The error of the cost function")
plt.show()
运行结果如下: