牛顿下降回溯直线及基于停止准则的搜索算法结果分析和比较（XJTU优化课第四章算法实现）

今天这个程序的主要是比较在求一个约束函数的最优值时，采取停止准则的搜索方法和回溯直线搜索方法的优劣。话不多说，先上结论：就我所测试的函数，我改变了多个测试函数和不同的初始值，并测试了不同的回溯因子。结论是我要吹爆停止准则的搜索方法，该方法总是能够获得较好的优化值，但是回溯法甚至不能满足初始迭代条件，即使能够满足，得到的优化结论误差也大得离谱。话不多说，我直接上相关算法，算法如下：

基于以上算法，我构造了一个二元函数，并编写了以下代码进行测试。读者可以直接上手我的代码进行测试。

#等式约束的牛顿方法+回溯直线搜索实现
#结束条件选择停止准则和回溯搜索条件准则，并将这两种计算的结果进行比较

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

def function(x1,x2):
    value=math.pow(x1,2)-x1*x2+0.5*math.pow(x2,2)-2*x1+math.exp(x1-x2)
    return value 
def Jac_x1(x1,x2):
    value=2*x1-x2+math.exp(x1-x2)
    return value
def Jac_x2(x1,x2):
    value=-x1+x2-math.exp(x1-x2)
    return value

#设置初始数据值
X=[1,-3] #设置初始值 X的选择应该满足等式约束条件
x=[X[0],X[1]] #记录初始点
Jac=[Jac_x1(x[0],x[1]),Jac_x2(x[0],x[1])] #求解一阶雅可比矩阵
Hessian=np.array([[2,-1],[-1,1]])
Hessian=np.matrix(Hessian)
Hessian_I=Hessian.I #Hessian矩阵求逆
Direction=[-(Hessian_I[0,0]*Jac[0]+Hessian_I[0,1]*Jac[1]),-(Hessian_I[1,0]*Jac[0]+Hessian_I[1,1]*Jac[1])] #计算牛顿方向初始值
Lambda=(Direction[0]*Hessian[0,0]+Direction[1]*Hessian[1,0])*Direction[0]+(Direction[0]*Hessian[0,1]+Direction[1]*Hessian[1,1])*Direction[1]  #计算 步长增量 停止准则的计算方程
Epsilon=0.000000001 #设置停止准则，即精度
print(Lambda)

#回溯搜索的参数设置
Alpha=0.01
Beta=0.8#设置回溯系数
Tk=1 #步长参数设置

#设置存储向量
result1=[]  #存储停止准则的结果
error1=[] #存储停止准则下的误差的向量
Minimum1=function(X[0],X[1]) #停止准则下最优化值初始化最小值
count1=0 #停止准则下初始化迭代次数


result2=[]  #存储回溯搜索下的结果
error2=[] #存储回溯搜索下的的误差的向量
Minimum2=function(X[0],X[1]) #回溯搜索下下最优化值初始化最小值
count2=0 #回溯搜索下下初始化迭代次数
while 0.5*Lambda>Epsilon: #t停止准则下的搜索过程
    #记录结果
    result1.append(function(X[0],X[1]))
    if(function(X[0],x[1])<Minimum1):
        Minimum1=function(X[0],X[1])
    count1+=1
    

    #更新相关点
    X[0]=X[0]+Tk*Direction[0]
    X[1]=X[1]+Tk*Direction[1]
    Jac=[Jac_x1(X[0],X[1]),Jac_x2(X[0],X[1])] #求解一阶雅可比矩阵
    Hessian=np.array([[2,-1],[-1,1]])
    Hessian=np.matrix(Hessian)
    Hessian_I=Hessian.I #Hessian矩阵求逆
    Direction=[-(Hessian_I[0,0]*Jac[0]+Hessian_I[0,1]*Jac[1]),-(Hessian_I[1,0]*Jac[0]+Hessian_I[1,1]*Jac[1])] #计算牛顿方向初始值
    Lambda=(Direction[0]*Hessian[0,0]+Direction[1]*Hessian[1,0])*Direction[0]+(Direction[0]*Hessian[0,1]+Direction[1]*Hessian[1,1])                 *Direction[1]  #计算 步长增量 停止准则的计算方程
for i in range(0,count1,1):
    error1.append(result1[i]-Minimum1)

#停止准则下的结果展示
print("初始迭代点是：({}，{})".format(x[0],x[1]))
print("停止准则下迭代的最优值是：{}".format(Minimum1))
print("停止准则下的迭代的最优点是：（{},{})".format(X[0],X[1]))
print("停止准则下迭代的次数是：{}".format(count1))
plt.plot(result1)
plt.title("The iteration process under the Limitation conditions")
plt.xlabel("The number of iteration")
plt.ylabel("The value of iteration")
plt.show()
plt.figure()
plt.plot(error1)
plt.title("The error under the limitation conditions")
plt.xlabel("The number of iteration")
plt.ylabel("The error in the process of iteration")
plt.show()

#恢复原始数据进行回溯方法测试
X=[x[0],x[1]] #恢复初始点准备使用回溯方法进行测试
Jac=[Jac_x1(x[0],x[1]),Jac_x2(x[0],x[1])] #求解一阶雅可比矩阵
Hessian=np.array([[2,-1],[-1,1]])
Hessian=np.matrix(Hessian)
Hessian_I=Hessian.I #Hessian矩阵求逆
Direction=[-(Hessian_I[0,0]*Jac[0]+Hessian_I[0,1]*Jac[1]),-(Hessian_I[1,0]*Jac[0]+Hessian_I[1,1]*Jac[1])] #计算牛顿方向初始值
Lambda=(Direction[0]*Hessian[0,0]+Direction[1]*Hessian[1,0])*Direction[0]+(Direction[0]*Hessian[0,1]+Direction[1]*Hessian[1,1])*Direction[1]  #计算 步长增量 停止准则的计算方程


print(function(X[0]+Tk*Direction[0],X[1]+Tk*Direction[1]))
print((function(X[0],X[1])-Alpha*Tk*0.5*Lambda))
print(function(X[0],X[1]))
print(Lambda)
while function(X[0]+Tk*Direction[0],X[1]+Tk*Direction[1])>(function(X[0],X[1])-Alpha*Tk*0.5*Lambda):#回溯方法指定了一个间隔，即新的搜索必须大于此间隔
    result2.append(function(X[0],X[1]))
    if(Minimum2>function(X[0],X[1])):
        Minimum2=function(X[0],X[1])
    count2+=1
    Tk=Beta*Tk #回溯点设置
    print(function(X[0],x[1]))
    #更新点
    Jac=[Jac_x1(x[0],x[1]),Jac_x2(x[0],x[1])] #求解一阶雅可比矩阵
    Hessian=np.array([[2,-1],[-1,1]])
    Hessian=np.matrix(Hessian)
    Hessian_I=Hessian.I #Hessian矩阵求逆
    Direction=[-(Hessian_I[0,0]*Jac[0]+Hessian_I[0,1]*Jac[1]),-(Hessian_I[1,0]*Jac[0]+Hessian_I[1,1]*Jac[1])] #计算牛顿方向初始值
    Lambda=(Direction[0]*Hessian[0,0]+Direction[1]*Hessian[1,0])*Direction[0]+(Direction[0]*Hessian[0,1]+Direction[1]*Hessian[1,1])                 *Direction[1]  #计算 步长增量 停止准则的计算方程
for i in range(0,count2,1):
    error2.append(result2[i]-Minimum2)
#回溯搜索下的迭代结果展示
plt.figure()
print("回溯搜索下迭代的最优值是：{}".format(Minimum2))
print("回溯搜索下的迭代的最优点是：（{},{})".format(X[0],X[1]))
print("回溯搜索下迭代的次数是：{}".format(count2))
plt.plot(result2)
plt.title("The iteration process under the traceback conditions")
plt.xlabel("The number of iteration")
plt.ylabel("The value of iteration")
plt.show()
plt.figure()
plt.plot(error2)
plt.title("The error under the traceback conditions")
plt.xlabel("The number of iteration")
plt.ylabel("The error in the process of iteration")
plt.show()

代码有点多哈哈，不过结构大概就是定义函数，赋初始值，进行两种搜索条件的迭代过程并记录结果，迭代过程伴随点的更新。最后将记录的结果展示出来，这就是这个程序的大致思路。以下展示这个程序运行的结果：

以上是采用停止准则获得的函数变化分析，可以看到它最终还是得到了一个比较好的优化值，迭代次数也可接受，它的误差函数图如下：

下面将展示回溯直线搜索的结果。
下图展示了回溯直线搜索的函数运行图和迭代过程：

好家伙，直接给爷一条直线，我cao了。好吧，哪位大哥看看可不可以改改对应参数使它搜索到更好的优化值，而不是60多这么离谱的值。停止准则的搜索结果大概是1，这个差距！下面我们来看看回溯直线搜索的误差图：

哈哈，不用想，这个回溯法让人大开眼界，这我直接日了，误差基本没啥变化，就这？？？好吧，也许爷没调好参，爷哭了，反正我是修改了很多次参数各种改，反正结果始终很离谱，不信老哥可以去试试。哈哈
最终一句话：回溯直线不太好调参。难受ing…