cover

最优化理论——牛顿法（一维、多维）原理及python代码实现

代码如下，如果只想抄个模板可以直接复制一维牛顿法代码二维牛顿法实现。

weixin_46668520

4254人浏览 · 2023-11-07 08:47:10

weixin_46668520 · 2023-11-07 08:47:10 发布

目录

1.一维牛顿法

1.1一维牛顿法原理

1.2一维牛顿法代码实现

2.多维牛顿法

2.1多维牛顿法原理

2.2多维牛顿法代码实现

代码如下，如果只想抄个模板可以直接复制

一维牛顿法代码

from  sympy import *
import pylab
x = symbols('x')
fun = eval(input("请输入一维函数"))
x_0 = float(input("请输入初始点"))
epsilon = float(input("输入精度"))

diff_1 = diff(fun)
diff_2 = diff(diff_1)

def newton(x_0,epsilon):
    while True:
        try:
             x_0 = x_0 - diff_1.subs(x,x_0) / diff_2.subs(x,x_0)
        except ZeroDivisionError:
            raise ZeroDivisionError("Could not find root") from None
        if abs(diff_1.subs(x,x_0)) < epsilon:
            return x_0
x_min = newton(x_0,epsilon)
print(x_min)
print("所求极值点为{:.2f}：".format(x_min))
print("所求极值为{:.5f}".format(fun.subs(x,x_min)))
def drawf(a,b,interp):
    x_la = [a + ele * interp for ele in range(0, int((b - a) / interp))]
    y_la = []
    for ele in x_la:
        y_la.append(fun.subs(x,ele).evalf())
    pylab.plt.figure(1)
    pylab.plt.plot(x_la, y_la)
    pylab.xlim(a, b)
    pylab.plt.show()
drawf(-10,10,0.05)

二维牛顿法实现

from  sympy import *
import pylab
import numpy as np
import math

x1 = symbols('x1')
x2 = symbols('x2')
fun = eval(input("请输入二维函数（用x1,x2定义）"))
x_0 = input("请输入初始点，用逗号隔开")
x_0 = x_0.split(',')
data = []
for i in x_0:
    data.append(float(i))
print(data)
epsilon = float(input("输入精度"))
def grad_len(grad):
    vec_len = math.sqrt(pow(grad[0], 2) + pow(grad[1], 2))
    return vec_len
def newton_dou(x_init, fun,epsilon):
    # 初始点的梯度值
    grandient_obj = np.array([diff(fun, x1).subs(x1, x_init[0]).subs(x2, x_init[1]),
                              diff(fun, x2).subs(x1, x_init[0]).subs(x2, x_init[1])], dtype=float)
    # 初始点的hessian矩阵
    hessian_obj = np.array([[diff(fun, x1, 2).subs(x1, x_init[0]).subs(x2, x_init[1]),
                             diff(diff(fun, x1), x2).subs(x1, x_init[0]).subs(x2, x_init[1])],
                            [diff(diff(fun, x2), x1).subs(x1, x_init[0]).subs(x2, x_init[1]),
                             diff(fun, x2, 2).subs(x1, x_init[0]).subs(x2, x_init[1])]],dtype=float)
    # 海森矩阵求逆
    inverse = np.linalg.inv(hessian_obj)
    #inverse = np.linalg.inv(hessian_obj)
    if(grad_len(grandient_obj)>epsilon):
        # 第一次迭代公式
        x_new = x_init - np.matmul(inverse, grandient_obj)
        #print(x_new)
            # print('迭代第%d次：%.5f' %(i, x_new))
        # 当前点的梯度值
        grandient_obj = np.array([diff(fun, x1).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                  diff(fun, x2).subs(x1, x_new[0]).subs(x2, x_new[1])], dtype=float)
        # 当前点的hessian矩阵
        hessian_obj = np.array([[diff(fun, x1, 2).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                 diff(diff(fun, x1), x2).subs(x1, x_new[0]).subs(x2, x_new[1])],
                                [diff(diff(fun, x2), x1).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                 diff(fun, x2, 2).subs(x1, x_new[0]).subs(x2, x_new[1])]],dtype=float)

        inverse = np.linalg.inv(hessian_obj)  # hessian矩阵求逆
        x_new = x_new - np.matmul(inverse, grandient_obj)  # 迭代公式
        #print(x_new)
    return x_new
x_min = newton_dou(data,fun,epsilon)
print("所求的最小值点:x1为{:.5f},x2为{:.5f}".format(x_min[0],x_min[1]))
print("所求的最小值为{:5f}".format(fun.subs({x1:x_min[0],x2:x_min[1]})))

1.一维牛顿法

一维牛顿法（One-Dimensional Newton's Method）是一种在最优化问题中使用的迭代算法，用于寻找函数的极值点。它是基于牛顿法的一维优化版本，适用于单变量函数。

1.1一维牛顿法原理

原理介绍：对于要求的函数q(x),我们使用泰勒展开式，假设q(x)在所求区间内二阶可导，可推导：

$q(x) = f(x_k) + f'(x_k)(x-x_k) + \frac{1}{2}f''(x)(x-x_{k})^2$

对q(x)求导得：

$q'(x) = f'(x_k) + f''(x_k)(x-x_0)$

要求极值，我们令q'(x) = 0，得x的表达式

$x = x_k - \frac{f'(x)}{f''(x)}$

由此，我们得到了迭代公式，我们只需要每次更新x，并且设置一个精度，当:

$|f'(x)|<\theta$

停止迭代.

1.2一维牛顿法代码实现

一维牛顿法代码实现比较简单，代码如下。

我们在此输入函数为5*x**2 +2*x +3 精度为0.001，初始点为10

$5x^2+2x +3$

from  sympy import *
import pylab
x = symbols('x')
fun = eval(input("请输入一维函数"))
x_0 = float(input("请输入初始点"))
epsilon = float(input("输入精度"))

diff_1 = diff(fun)
diff_2 = diff(diff_1)

def newton(x_0,epsilon):
    while True:
        try:
             x_0 = x_0 - diff_1.subs(x,x_0) / diff_2.subs(x,x_0)
        except ZeroDivisionError:
            raise ZeroDivisionError("Could not find root") from None
        if abs(diff_1.subs(x,x_0)) < epsilon:
            return x_0
x_min = newton(x_0,epsilon)
print(x_min)
print("所求极值点为{:.2f}：".format(x_min))
print("所求极值为{:.5f}".format(fun.subs(x,x_min)))

运行结果为

绘制原函数图像,这里我们定义了a，b和interp的取值分别为-10，10，0.05。大家可以根据自己需求进行修改。

def drawf(a,b,interp):
    x_la = [a + ele * interp for ele in range(0, int((b - a) / interp))]
    y_la = []
    for ele in x_la:
        y_la.append(fun.subs(x,ele).evalf())
    pylab.plt.figure(1)
    pylab.plt.plot(x_la, y_la)
    pylab.xlim(a, b)
    pylab.plt.show()
drawf(-10,10,0.05)

绘制的图像为

$x = x_k - \Delta ^2f(x_k)^{-1}\Delta f(x_k)$

2.多维牛顿法

2.1多维牛顿法原理

多维牛顿法，在形式上与一维牛顿法相似：

由泰勒公式求f(x)在x_k处的展开:

$f(x) = f(x_k) + \Delta ^tf(x_k)(x-x_k) + \frac{1}{2}(x-x_k)^t\Delta ^2f(x)(x-x_{k})$

求导得：要求海森矩阵为正定，可逆的，并且为对称矩阵

$\Delta f(x) + \Delta ^2f(x_k)(x-x_k) = 0$

得到迭代公式

$x = x_k - \Delta^2f(x)^{-1}\Delta f(x_k)$

算法设计流程

(1). 给定初始点 $x_1\in R^N$ 允许误差 $\theta$ >0

(2) 计算梯度，Hesse矩阵和牛顿方向

$\Delta f(x)$ $\Delta ^2f(x)$ $d_k = - \Delta^2f(x)^{-1}\Delta f(x_k)$

(3) 若 $||\Delta f(x)|| < \theta$ 停止计算，输出 $x^* = x$

否则：牛顿迭代

$x = x_k - \Delta^2f(x)^{-1}\Delta f(x_k)$

牛顿法虽然简单，收敛速度快（能够达到二次收敛），但也有许多缺点

（1）局部收敛性，收敛点必须非常接近最优解

（2）计算复杂度高

（3）hesse矩阵可能为奇异矩阵，其逆矩阵不一定存在

2.2多维牛顿法代码实现

导入相关的包

from  sympy import *
import pylab
import numpy as np
import math

输入和定义相关的参数，或者函数

输入的函数为 x1**2 + x1*x2*5 +x2**4

$x_1^2 +5x_1x_2+x_2^4$

初始点为 10，10 （注意逗号必须是英文符号）

x1 = symbols('x1')
x2 = symbols('x2')
fun = eval(input("请输入二维函数（用x1,x2定义）"))
x_0 = input("请输入初始点，用逗号隔开")
x_0 = x_0.split(',')
data = []
for i in x_0:
    data.append(float(i))
print(data)
epsilon = float(input("输入精度"))

精度设为0.0001

定义模函数，这里是二范数

def grad_len(grad):
    vec_len = math.sqrt(pow(grad[0], 2) + pow(grad[1], 2))
    return vec_len

定义牛顿方法

def newton_dou(x_init, fun,epsilon):
    # 初始点的梯度值
    grandient_obj = np.array([diff(fun, x1).subs(x1, x_init[0]).subs(x2, x_init[1]),
                              diff(fun, x2).subs(x1, x_init[0]).subs(x2, x_init[1])], dtype=float)
    # 初始点的hessian矩阵
    hessian_obj = np.array([[diff(fun, x1, 2).subs(x1, x_init[0]).subs(x2, x_init[1]),
                             diff(diff(fun, x1), x2).subs(x1, x_init[0]).subs(x2, x_init[1])],
                            [diff(diff(fun, x2), x1).subs(x1, x_init[0]).subs(x2, x_init[1]),
                             diff(fun, x2, 2).subs(x1, x_init[0]).subs(x2, x_init[1])]],dtype=float)
    # 海森矩阵求逆
    inverse = np.linalg.inv(hessian_obj)
    #inverse = np.linalg.inv(hessian_obj)
    if(grad_len(grandient_obj)>epsilon):
        # 第一次迭代公式
        x_new = x_init - np.matmul(inverse, grandient_obj)
        #print(x_new)
            # print('迭代第%d次：%.5f' %(i, x_new))
        # 当前点的梯度值
        grandient_obj = np.array([diff(fun, x1).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                  diff(fun, x2).subs(x1, x_new[0]).subs(x2, x_new[1])], dtype=float)
        # 当前点的hessian矩阵
        hessian_obj = np.array([[diff(fun, x1, 2).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                 diff(diff(fun, x1), x2).subs(x1, x_new[0]).subs(x2, x_new[1])],
                                [diff(diff(fun, x2), x1).subs(x1, x_new[0]).subs(x2, x_new[1]),
                                 diff(fun, x2, 2).subs(x1, x_new[0]).subs(x2, x_new[1])]],dtype=float)

        inverse = np.linalg.inv(hessian_obj)  # hessian矩阵求逆
        x_new = x_new - np.matmul(inverse, grandient_obj)  # 迭代公式
        #print(x_new)
    return x_new

调用相关函数

x_min = newton_dou(data,fun,epsilon)
print("所求的最小值点:x1为{:.5f},x2为{:.5f}".format(x_min[0],x_min[1]))
print("所求的最小值为{:5f}".format(fun.subs({x1:x_min[0],x2:x_min[1]})))

得到运行结果

腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

cover

动物识别系统 python实现+UI

腾讯云开发者社区

Linux下安装MySQL8.0(超详细)

Linux下安装MySQL8.0(超详细)

腾讯云开发者社区

cover

安卓手机如何连接苹果电脑

腾讯云开发者社区

所有评论(0)

查看更多评论

weixin_46668520

@weixin_46668520

已为社区贡献2条内容