最优化理论 无约束规划 Newton法 牛顿法 求极小点

用到 numpy和sympy 库

问题

对于目标函数如下，用 Newton 法求其极小点。代码中我设置的精度为10^(-3)

Newton法基本步骤

①取初始点x(1) ，初始化精度e，设置迭代次数 k = 1；
②当递归边界 || ▽f(x(k)) || <= e（其实就是梯度值差不多接近于0的时候，趋近于极小值，类比单变量的二次函数对 x* 求导为0的地方为极值），当前x(k)作为最优解；否则，解方程得到 d(k)；

③ 置 x(k+1) = x(k) + d ， k += 1，转步骤②。
牛顿法的相关证明晚些我会整理出来写一篇文章。

代码

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
# f = x1 ** 2 + x2 ** 2 + x1 * x2 + 2 * x1 - 3 * x2
f = 4 * x1**2 + x2**2 - (x1**2) * x2  # 目标函数
x0 = np.array([1, 1])  # 初始点
e = 1e-3  # 精度


def Newton(x, k):
    df1 = sp.diff(f, x1)  # 对x1求偏导
    df2 = sp.diff(f, x2)  # 对x2求偏导
    # g_df = np.array([[df1], [df2]])  # 获取一阶梯度向量
    # 将值代入一阶梯度中
    g_dfx = np.array([df1.evalf(subs={x1: x[0], x2: x[1]}), df2.evalf(subs={x1: x[0], x2: x[1]})], dtype=np.float64)  # 获得已被赋值的梯度
    # print(g_dfx)
    d2f1 = sp.diff(df1, x1)  # 对x1求二阶偏导，下同
    d2f2 = sp.diff(df1, x2)
    d2f1_ = sp.diff(df2, x1)
    d2f2_ = sp.diff(df2, x2)
    # g_d2f = np.array([[d2f1, d2f2], [d2f1_, d2f2_]])
    # 获取代入值的汉森矩阵
    g_d2fx = np.array([[d2f1.evalf(subs={x1: x[0], x2: x[1]}),
                        d2f2.evalf(subs={x1: x[0], x2: x[1]})],
                       [d2f1_.evalf(subs={x1: x[0], x2: x[1]}),
                        d2f2_.evalf(subs={x1: x[0], x2: x[1]})]], dtype=np.float64)
    # print(g_d2fx)
    dk = np.dot(-g_dfx, np.linalg.inv(g_d2fx))
    xk = x + dk
    k += 1
    print("(%d)  x =" % k, x, "dk =", dk)
    print(g_dfx)
    print(g_d2fx)
    if np.linalg.norm(g_dfx) <= e:
        print("范数值为：", np.linalg.norm(g_dfx))
        print("最终解为：", x)
        return
    else:
        Newton(xk, k)


Newton(x0, 0)

结尾

其实开头那个目标函数迭代次数很少，会发现迭代一次就结束了，所以可以试试代码中被注释的那个函数，迭代次数会多一些，能更直观感受牛顿法求解极小点的过程。