编程实现newton插值c++_数值最优化方法-Newton方法

参考文献：高立. 数值最优化方法[M]. 北京大学出版社, 2014.

本文主要基于Python实现第三章的Newton算法，对于理论部分的推导详见参考文献。

考虑问题

从不同的初始点出发用Newton方法求解该问题。起始点的坐标值分别为

，

，

。

先用Python画出等高线图和3D图看一下函数的样子，极小点应该大致位于下图五角星的位置。

下面使用基本Newton法和LM算法来求解该问题。

算法1.基本Newton方法

基本Newton方法指取全步长

的Newton方法。

基本Newton方法的迭代步骤：
1.给出

,

,

2.若终止准则满足，则输出有关信息，停止迭代
3.由方程组

计算

，式中

和

分别为

在第k次迭代时二阶导数和一阶导数的值。

（注：此时需满足

正定，划重点！！！）

4.

,

,转步2

Python代码实现如下，实现时进行了可视化来显示每次迭代过程，并可以输出第k次的迭代点坐标值以及函数值。

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


def get_Hessian(data):
    dx1x1 = y.diff(x1,x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx1x2 = y.diff(x1,x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2x1 = y.diff(x2,x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2x2 = y.diff(x2,x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    return np.array([[dx1x1, dx1x2],[dx2x1, dx2x2]])


def get_Jacobian(data):
    dx1 = y.diff(x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2 = y.diff(x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    return np.array([[dx1],[dx2]]), math.sqrt(dx1**2+dx2**2)


def plot_contour():
    x = np.linspace(-6,6,100)
    y = np.linspace(-6,6,100)
    X,Y = np.meshgrid(x,y)
    Z = 3*X**2+3*Y**2-X**2*Y
    plt.contour(X, Y, Z, 25, colors = 'red', linewidth = 0.5)


if __name__ == "__main__":
    # 打开交互模式
    plt.ion()
    # 绘制等高线
    plot_contour()
    plt.xlabel('x1')
    plt.ylabel('x2')
    # 定义迭代过程需求解的函数
    x1,x2,y = symbols("x1 x2 y")
    y = 3*x1**2+3*x2**2-x1**2*x2
    # 定义起始点
    x_initial = [1.5,1.5]
    x_k=x_initial
    xx = []
    yy = []
    zz = []
    xx.append(round(x_k[0],4))
    yy.append(round(x_k[1],4))
    zz.append(round(y.subs(x1,x_k[0]).subs(x2,x_k[1]).evalf(),4))
    print("({:.4f},{:.4f}),    {:.4f}".format(xx[0], yy[0], zz[0]))
    plt.plot(xx,yy,c='black')
    plt.scatter(xx,yy,s=30,c='b')
    plt.pause(5)
    for i in range(1,7):
        Gk = get_Hessian(x_k)
        Gk = Gk.astype(float)
        gk, gknorm = get_Jacobian(x_k)
        gk = gk.astype(float)
        Gk_inverse = np.linalg.inv(Gk)
        dk = np.dot(Gk_inverse, gk)
        dk = dk.T[0]
        x_k = [x_k[0]-dk[0],x_k[1]-dk[1]]
        xx.append(round(x_k[0],4))
        yy.append(round(x_k[1],4))
        zz.append(round(y.subs(x1,x_k[0]).subs(x2,x_k[1]).evalf(),4))
        print("({:.4f},{:.4f}),    {:.4f}".format(xx[i], yy[i], zz[i]))
        plt.plot(xx,yy,c='black')
        plt.scatter(xx,yy,s=30,c='b')
        plt.pause(5)
        if gknorm <= 10e-6:
            break
    plt.ioff()
    plt.show()   

对于起始点为

，迭代过程如下，经过千辛万苦最终还是到达了五角星的附近。

图中，红色线为等高线

迭代点及其函数值

起始点为

，其最终迭代结果如下，这次并没有太幸运，没有胜利到达五角星。

起始点为

时，由于此时

为奇异阵，此时基本Newton法失败。干脆就直接罢工了。由此可见，基础Newton法求最优化问题对于初始点的选择要求较高，初始点选取的不好，算法结果会较差。

在迭代过程中出现G非正定的情况亦会导致Newton法失败。算法2的LM方法可以解决该问题。

算法2.LM方法

LM方法即Levenberg-Marquardt方法，是处理

奇异、不正定等情形的一个最简单且有效的方法。该方法将基本Newton法第3步求解步长的方程

更改为

式中

,

为单位阵。当

也不正定时，简单令

即可。

Python实现如下：

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


def get_Hessian(data):
    dx1x1 = y.diff(x1,x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx1x2 = y.diff(x1,x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2x1 = y.diff(x2,x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2x2 = y.diff(x2,x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    return np.array([[dx1x1, dx1x2],[dx2x1, dx2x2]])


def get_Jacobian(data):
    dx1 = y.diff(x1).subs(x1,data[0]).subs(x2,data[1]).evalf()
    dx2 = y.diff(x2).subs(x1,data[0]).subs(x2,data[1]).evalf()
    return np.array([[dx1],[dx2]]), math.sqrt(dx1**2+dx2**2)


def plot_contour():
    x = np.linspace(-6,6,100)
    y = np.linspace(-6,6,100)
    X,Y = np.meshgrid(x,y)
    Z = 3*X**2+3*Y**2-X**2*Y
    plt.contour(X, Y, Z, 25, colors = 'red', linewidth = 0.5)


def judge_positive(G_k):
    # 计算特征值，当特征值均大于0时矩阵正定
    eigenvalue,featurevector=np.linalg.eig(G_k)
    for i in eigenvalue:
        if i <= 0:
            return False
    return True


if __name__ == "__main__":
    # 打开交互模式
    plt.ion()
    # 绘制等高线
    plot_contour()
    plt.xlabel('x1')
    plt.ylabel('x2')
    # 定义迭代过程需求解的函数
    x1,x2,y = symbols("x1 x2 y")
    y = 3*x1**2+3*x2**2-x1**2*x2
    # 定义起始点
    x_initial = [-2,4]
    x_k=x_initial
    # 定义vk
    vk = 3.02
    xx = []
    yy = []
    zz = []
    xx.append(round(x_k[0],4))
    yy.append(round(x_k[1],4))
    zz.append(round(y.subs(x1,x_k[0]).subs(x2,x_k[1]).evalf(),4))
    print("({:.4f},{:.4f}),    {:.4f}".format(xx[0], yy[0], zz[0]))
    plt.plot(xx,yy,c='black')
    plt.scatter(xx,yy,s=30,c='b')
    plt.pause(1)
    for i in range(1,7):
        Gk = get_Hessian(x_k)
        Gk = Gk.astype(float)
        I = np.identity(2)
        Gk_plus = Gk + vk*I
        # 判断Gk + vk*I的正定性
        if judge_positive(Gk_plus):
            Gk = Gk_plus
        else:
            Gk = Gk + 2*vk*I
        gk, gknorm = get_Jacobian(x_k)
        gk = gk.astype(float)
        Gk_inverse = np.linalg.inv(Gk)
        dk = np.dot(Gk_inverse, gk)
        dk = dk.T[0]
        x_k = [x_k[0]-dk[0],x_k[1]-dk[1]]
        xx.append(round(x_k[0],4))
        yy.append(round(x_k[1],4))
        zz.append(round(y.subs(x1,x_k[0]).subs(x2,x_k[1]).evalf(),4))
        print("({:.4f},{:.4f}),    {:.4f}".format(xx[i], yy[i], zz[i]))
        plt.plot(xx,yy,c='black')
        plt.scatter(xx,yy,s=30,c='b')
        plt.pause(1)
        if gknorm <= 10e-6:
            break
    plt.ioff()
    plt.show()

在Python的实现中，对于矩阵

的判定采用特征值法，即求出矩阵

的特征值，当该矩阵的特征值均为正时，矩阵正定。

LM算法的迭代过程如下图所示，

取3.02。可以看出当起始点为

时也能正常迭代。起始点为

和

时，迭代效果要好于基本Newton法。三次结果均成功抵达了五角星。