python实现最速下降法和牛顿法

小白上路，将每次学习到的一点点积累在这，便于以后查阅。

问题：

使用回溯线性搜索方法实现最速下降法和牛顿法的编程，其中初始点：x0=(-1.2,1)，初始步长：alpha0=1

预备知识：

（1） 回溯线性搜索

（2）最速下降法和牛顿法

（这是我的笔记有点潦草，懒得再写，就直接贴图过来了，重点是绿色标记的部分）

编程代码：

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.axisartist as axisartist

#回溯线性搜索   
def BLS(f,df,xk,alpha,pk):
    t=0.5
    c=1e-4
    f_k=f.subs(zip(x,xk))
    df_k=[df[0].subs(zip(x,xk)),df[1].subs(zip(x,xk))]
    while f.subs(zip(x,xk+alpha*pk))>f_k+c*alpha*np.dot(df_k,pk):
        alpha=alpha*t
    return alpha
 
 #最速下降法   
def steepest_descent(x0,alpha0,f,df):
    k=0
    xk=x0
    data_x1=[x0[0]]
    data_x2=[x0[1]]
    alpha=alpha0
    alpha_save=[alpha]
    df_k=[df[0].subs(zip(x,xk)),df[1].subs(zip(x,xk))]
    df_k=np.array(df_k,dtype=np.float64)
    while np.linalg.norm(df_k)>1e-6:
        pk=[-1/np.linalg.norm(df_k)]*df_k
        alpha=BLS(f,df,xk,alpha,pk)
        alpha_save.append(alpha)
        xk=xk+alpha*pk
        k=k+1
        data_x1.append(xk[0])
        data_x2.append(xk[1])
        df_k=[df[0].subs(zip(x,xk)),df[1].subs(zip(x,xk))]
        df_k=np.array(df_k,dtype=np.float64)
    return data_x1,data_x2,k,alpha_save
    
#牛顿法
def Newton(x0,alpha0,f,df,hess):
    k=0
    xk=x0
    data_x1=[x0[0]]
    data_x2=[x0[1]]
    alpha=alpha0
    alpha_save=[alpha]
    df_k=[df[0].subs(zip(x,xk)),df[1].subs(zip(x,xk))]
    df_k=np.array(df_k,dtype=np.float64)
    hess_k=[[hess[0,0].subs(zip(x,xk)),hess[0,1].subs(zip(x,xk))],[hess[1,0].subs(zip(x,xk)),hess[1,1].subs(zip(x,xk))]]
    hess_k=np.array(hess_k,dtype=np.float64)
    while np.linalg.norm(df_k)>1e-6:
        pk=np.dot(-1*np.linalg.inv(hess_k),df_k)
        alpha=BLS(f,df,xk,alpha,pk)
        alpha_save.append(alpha)
        xk=xk+alpha*pk
        k=k+1
        data_x1.append(xk[0])
        data_x2.append(xk[1])
        df_k=[df[0].subs(zip(x,xk)),df[1].subs(zip(x,xk))]
        df_k=np.array(df_k,dtype=np.float64)
        hess_k=[[hess[0,0].subs(zip(x,xk)),hess[0,1].subs(zip(x,xk))],[hess[1,0].subs(zip(x,xk)),hess[1,1].subs(zip(x,xk))]]
        hess_k=np.array(hess_k,dtype=np.float64)
    return data_x1,data_x2,k,alpha_save 
         
if __name__=="__main__":
    #初始化
    x1=symbols('x1')
    x2=symbols('x2')
    x=[x1,x2]
    x0=[-1.2,1]
    alpha0=1

    #求解最终解和迭代次数
    f=100*(x2-x1**2)**2+(1-x1)**2
    df=[diff(f,x1),diff(f,x2)]
    hess=np.array([[diff(df[0],x1),diff(df[0],x2)],[diff(df[1],x1),diff(df[1],x2)]])
    data_x1_s,data_x2_s,k_s,alpha_s=steepest_descent(x0,alpha0,f,df)
    X_s=[data_x1_s[-1],data_x2_s[-1]]
    print("BLS-steepest descent 最终解：",X_s)
    print("迭代次数：",k_s)
    data_x1_n,data_x2_n,k_n,alpha_n=Newton(x0,alpha0,f,df,hess)
    X_n=[data_x1_n[-1],data_x2_n[-1]]
    print("BLS-Newton 最终解：",X_n)
    print("迭代次数：",k_n)

    #绘图：迭代过程
    fig = plt.figure(1)
    ax = axisartist.Subplot(fig, 111)
    fig.add_axes(ax)
    ax.axis["bottom"].set_axisline_style("-|>", size=1.5)
    ax.axis["left"].set_axisline_style("->", size=1.5)
    ax.axis["top"].set_visible(False)
    ax.axis["right"].set_visible(False)
    plt.title(r'$BLS \ -\ steepest \ descent \ method$')
    plt.plot(data_x1_s, data_x2_s, label=r'$f=100 \cdot(x2-x1^2)^2+(1-x1)^2$')
    plt.legend()
    plt.scatter(data_x1_s[0], data_x2_s[0], s=100, c=20, marker=(9,2, 30), )
    plt.scatter(1, 1, marker=(5, 1), c=5, s=1000)
    plt.grid()
    plt.xlabel(r'$x_1$', fontsize=20)
    plt.ylabel(r'$x_2$', fontsize=20)
    plt.show()

    fig = plt.figure(2)
    ax = axisartist.Subplot(fig, 111)
    fig.add_axes(ax)
    ax.axis["bottom"].set_axisline_style("-|>", size=1.5)
    ax.axis["left"].set_axisline_style("->", size=1.5)
    ax.axis["top"].set_visible(False)
    ax.axis["right"].set_visible(False)
    plt.title(r'$BLS \ -\ Newton \ method$')
    plt.plot(data_x1_n, data_x2_n, label=r'$f=100 \cdot(x2-x1^2)^2+(1-x1)^2$')
    plt.legend()
    plt.scatter(data_x1_n[0], data_x2_n[0], s=100, c=20, marker=(9,2, 30), )
    plt.scatter(1, 1, marker=(5, 1), c=5, s=1000)
    plt.grid()
    plt.xlabel(r'$x_1$', fontsize=20)
    plt.ylabel(r'$x_2$', fontsize=20)
    plt.show()

    #绘图：迭代误差
    fig = plt.figure(3)
    xx=np.array([1,1])
    ks=[]
    errs=[]
    kn=[]
    errn=[]
    for i in range(len(data_x1_s)):
        ks.append(i)
        errs.append(np.linalg.norm(np.array([data_x1_s[i],data_x2_s[i]])-xx))
    for i in range(len(data_x1_n)):
        kn.append(i)
        errn.append(np.linalg.norm(np.array([data_x1_n[i],data_x2_n[i]])-xx))
    plt.title(r'$error$')
    plt.plot(ks, errs, color='blue',label=r'$steepest \ descent$')
    plt.plot(kn, errn, color='red',label=r'$Newton $')
    plt.yscale('log')
    plt.legend()
    plt.grid()
    plt.xlabel(r'$iterations$', fontsize=20)
    plt.ylabel(r'$error$', fontsize=20)
    plt.ylim(1e-7,1e+3)
    plt.show()

运行结果：

（每一次迭代后的xk这里没有写进来，实在太多了）

最速下降法：
最终解：[0.9999993331952023,
0.9999986637225776]
迭代次数：15671

牛顿法:
最终解：[0.9999998896294611,
0.9999997768698554]
迭代次数：174