一元和多元函数求极值(Python)-牛顿法

作者: 一个扫地的垃圾 | 来源:发表于2019-04-13 14:18 被阅读0次

一元和多元函数求极值(Python)-牛顿法
无约束凸优化算法
牛顿法求极值
拟牛顿法的原理
拉格朗日乘子法
无约束非线性规划的极值条件
多元函数的极值及其求法
Levenberg-Maquardt Algorithm 推导
支持向量机
极值理论

1. 引言

我们在中学的时候学过一元二次函数，求解时引入一个求根公式，代入公式就可以得到不同的根，假如想计算一个高次方程的解，我们还能推导出求根公式吗？
伽罗瓦在群论中证实，五次及以上多项式方程没有显示表达的求根公式，但是科学研究（例如行星轨道计算）中还是有很多求解高次方程的真实需求。既然得不到明确的求根公式，我们可以用迭代的办法来不断逼近真实的解。
多项式方程求解的问题实际上可以看成是函数求极值，假如我们对 $f'(x)=0$ 的求解得到驻点，将驻点代入 $f(x)$ 就可以计算出函数的极值。显然低次方程的求解可以用求根公式精确计算 $x$ 的值，但是高次方程的求解就要用逼近的思路比如梯度下降法、最速下降法、牛顿法等。本文主要讨论牛顿法和使用python对一元函数和多元函数求极值。

2. 牛顿法基本原理

2.1. 一元函数牛顿法

假如函数 $f(x)$ 是一个一元函数，使用泰勒公式二阶近似可以得到

$f(x)\approx f(x_0)+f'(x_0)(x-x_0)+\frac{1}{2}f''(x_0)(x-x_0)^2$

一元函数取极值的条件是 $f'(x)=0$ ，两端分别求导（下式约等号写作等号）得到

$f'(x)=f(x_0)+f''(x_0)(x-x_0)=0$

整理出 $x$ 的迭代公式

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

进一步写出迭代关系表达式，得到

$x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}$

因为牛顿法的迭代关系是根据要求解极值函数的导数图像在该点的斜率进行迭代，因此没有没有步长的概念，在哪一点得到极值仅仅依靠迭代的步数。相应的，迭代的次数越多越接近极值点。
上述在介绍了抽象的公式之后，我们来看这样一个小例子，已知一个函数 $f(x)=x^6+x$ ，在没有其他约束条件求 $f(x)$ 的极值。
在用牛顿法计算极值之前，使用软件绘制 $f(x)$ 的图像，用肉眼粗略的观测一下极值在哪里，可以帮助我们在后续计算中检验正确与否。

image.png

首先，需要给一个初始值 $x_0=10$ ,这是要开始进行迭代的位置；
然后，计算初始位置的一阶导数和二阶导数值， $f'(10)$ 和 $f''(10)$ ；
其次，代入迭代公式 $x_{1}=x_0-\frac{f'(x_0)}{f''(x_0)}$ ;
最后，将得到的第一次迭代的值 $x_1$ 重复上述两个步骤，得到 $x_2$ 、 $x_3$ 、··· $x_n$ ，停止条件为所设定的迭代步数。

2.2. 一元函数牛顿法的python程序

from sympy import *
# step为迭代步数，x0为初始位置，obj为要求极值的函数
def newtons(step, x0, obj):
    i = 1 # 记录迭代次数的变量
    x0 = float(x0) # 浮点数计算更快
    obj_deri = diff(obj, x) # 定义一阶导数，对应上述公式
    obj_sec_deri = diff(obj, x, 2) # 定义二阶导数，对应上述公式
    while i <= step:
        if i == 1:
            # 第一次迭代的更新公式
            xnew = x0 - (obj_deri.subs(x, x0)/obj_sec_deri.subs(x, x0))
            print('迭代第%d次：%.5f' %(i, xnew))
            i = i + 1
        else:
            #后续迭代的更新公式
            xnew = xnew - (obj_deri.subs(x, xnew)/obj_sec_deri.subs(x, xnew))
            print('迭代第%d次：%.5f' % (i, xnew))
            i = i + 1
    return xnew
x = symbols("x") # x为字符变量
result = newtons(50, 10, x**6+x)
print('最佳迭代的位置：%.5f' %result)

运行程序结果如下：

迭代第1次：8.00000
迭代第2次：6.39999
迭代第3次：5.11997
迭代第4次：4.09593
迭代第5次：3.27662
迭代第6次：2.62101
迭代第7次：2.09610
迭代第8次：1.67515
迭代第9次：1.33589
迭代第10次：1.05825
迭代第11次：0.82002
迭代第12次：0.58229
迭代第13次：0.17590
迭代第14次：-34.68063
迭代第15次：-27.74450
迭代第16次：-22.19560
迭代第17次：-17.75648
迭代第18次：-14.20519
迭代第19次：-11.36415
迭代第20次：-9.09132
迭代第21次：-7.27306
迭代第22次：-5.81846
迭代第23次：-4.65480
迭代第24次：-3.72391
迭代第25次：-2.97930
迭代第26次：-2.38386
迭代第27次：-1.90812
迭代第28次：-1.52901
迭代第29次：-1.22931
迭代第30次：-0.99804
迭代第31次：-0.83203
迭代第32次：-0.73518
迭代第33次：-0.70225
迭代第34次：-0.69886
迭代第35次：-0.69883
迭代第36次：-0.69883
迭代第37次：-0.69883
迭代第38次：-0.69883
迭代第39次：-0.69883
迭代第40次：-0.69883
迭代第41次：-0.69883
迭代第42次：-0.69883
迭代第43次：-0.69883
迭代第44次：-0.69883
迭代第45次：-0.69883
迭代第46次：-0.69883
迭代第47次：-0.69883
迭代第48次：-0.69883
迭代第49次：-0.69883
迭代第50次：-0.69883
最佳迭代的位置：-0.69883

Process finished with exit code 0

函数的极值点为-0.69883，与我们之前绘制的图像观测值一致。

2.3. 多元函数牛顿法

多元函数用牛顿法求极值的思路依然采用泰勒公式近似展开，只不过这里的 $X$ 为向量形式，假设 $f(X)$ 为二元函数:

$f(X)\approx f(X^{(0)}) +\nabla f(X^{(0)})^T(X-X^{(0)}) + \frac{1}{2}(X-X^{(0)})^T\nabla^2f(X^{(0)})(X-X^{(0)})$

对上述公式两侧求梯度，令其等于零，之后再进行整理，可得迭代公式(由于涉及矩阵求导数等一系列复杂操作，此处可以粗浅的理解一下)：

$X= X^{(0)} -H^{-1}f(X^{(0)}) \nabla f(X^{(0)})$

上述第一个公式中 $\nabla f(X^{(0)})$ 表示在 $X^{(0)}$ 点处的梯度值， $\nabla^2f(X^{(0)})$ 表示该点的二阶梯度也可以写成Hessian矩阵，迭代公式使用了Hessian矩阵的逆和梯度值来确定。

2.3. 多元函数牛顿法的python程序

from sympy import *
import numpy as np
# 假设多元函数是二维形式
# x_init为二维向量(x1, x2)
def newton_dou(step, x_init, obj):
    i = 1 # 记录迭代次数的变量
    while i <= step:
        if i == 1:
            grandient_obj = np.array([diff(obj, x1).subs(x1, x_init[0]).subs(x2, x_init[1]), diff(obj, x2).subs(x1, x_init[0]).subs(x2, x_init[1])], dtype=float) # 初始点的梯度值
            hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 初始点的hessian矩阵
            inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆
            x_new = x_init - np.matmul(inverse, grandient_obj) # 第一次迭代公式
            print(x_new)
            # print('迭代第%d次：%.5f' %(i, x_new))
            i = i + 1
        else:
            grandient_obj = np.array([diff(obj, x1).subs(x1, x_new[0]).subs(x2, x_new[1]), diff(obj, x2).subs(x1, x_new[0]).subs(x2, x_new[1])], dtype=float) # 当前点的梯度值
            hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 当前点的hessian矩阵
            inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆
            x_new = x_new - np.matmul(inverse, grandient_obj) # 迭代公式
            print(x_new)
            # print('迭代第%d次：%.5f' % (i, x_new))
            i = i + 1
    return x_new

x0 = np.array([0, 0], dtype=float)
x1 = symbols("x1")
x2 = symbols("x2")
newton_dou(5, x0, x1**2+2*x2**2-2*x1*x2-2*x2)

程序执行的结果如下：

[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]

Process finished with exit code 0

经过实际计算函数 $f(x)=x_1^2+2x_2^2-2x_1x_2-2x_2$ 的极值点为 $(1,1)$ ，只需一次迭代就能收敛到极值点。

3. 结束语

本文只是粗浅的将牛顿法的原理进行语言和程序描述，没有过多的探讨具体实施的细节，实际中机器学习等相关算法进行牛顿法迭代时并非如此编程。牛顿法由牛顿本人发现并命名，但是拉弗森早于牛顿发现的46年就已经提出该算法，现在牛顿法又称为牛顿-拉弗森方法。牛顿法有很多的缺点在此并没有讨论，例如在一元函数情况时，会产生导数为零不能迭代的情况或者迭代点有震荡往复等情况出现；多元函数的hessian矩阵计算困难，收敛过慢等情况，因此又产生了拟牛顿法和阻尼牛顿法等改进的方法，想进一步学习的读者可以参考非线性优化领域的相关书籍。

一元和多元函数求极值(Python)-牛顿法
1. 引言我们在中学的时候学过一元二次函数，求解时引入一个求根公式，代入公式就可以得到不同的根，假如想计...
无约束凸优化算法
本章涉及知识点1、scipy库求解全局最优和局最优2、多元函数的极值求解算法3、牛顿迭代法算法4、牛顿迭代法求解多...
牛顿法求极值
摘要最近在做两幅图像间单应性矩阵估计的研究，遇到LM（Levenberg-Marquardt）算法求极值，在查找...
拟牛顿法的原理
多元函数的泰勒展开式image-20200403212859301.png 牛顿法牛顿法是梯度下降法的进一步发展，...
拉格朗日乘子法
问题：求函数在条件下可能的极值点，其中 . 利用 Lagrange 乘子法，可将带约束的极值问题转化为无约束...
无约束非线性规划的极值条件
（一元函数极值点存在的必要条件和充分条件）对于二阶可微的一元函数，f(x)极值存在的必要条件是f'(x)=0 充分...
多元函数的极值及其求法
一元函数极值的必要条件二元函数极值的必要条件二元函数级值的充分条件 n 元函数极值的充分条件
Levenberg-Maquardt Algorithm 推导
前置知识 1. 牛顿法作用：1. 求根 2.求极值求根目标: 求解的根计算穿过初始点并且斜率为的直线...
支持向量机
1.问题求解： (1)拉格朗日乘子法定义拉格朗日函数，KKT条件为：求极值，则令得到：代入消去和,得到原问题的对...
极值理论
一元极值理论经典极值理论极值分布的类型以及性质极值分布的最大值吸引场平均超出量函数与T年重现水平平均超出...