最优化技术——牛顿法

发表于 | 分类于 |

最优化技术——牛顿法

没错，到期末了，我又开始学最优化了

参考目录

解决的问题

解决的是无约束的极小化问题，用数学语言描述如下： \[ \mathop{min}\limits_{x}f(x) \tag{0} \]

简单一维牛顿法

从现有极小值估计点出发，对\(f(x)\)做泰勒展开，进而找到极小值的下一个估计值。设\(x_k\)是当前极小值的估计值，有： \[ \varphi (x) = f(x_k) +f'(x_k)(x - x_k) + \frac{1}{2}f''(x_k)(x-x_k)^2 \tag{1} \] \(\varphi\)表示的是\(f(x)\) 在\(x\)附近的二阶泰勒展开式，由于求的是最值由极值条件可知:\(\varphi\)应当满足: \[ \varphi'(x) = 0 \tag{2} \] 将1、2两个式子整合可知： \[ f'(x_k) + f''(x_k)(x - x_k) = 0 \] 所以可以求得： \[ x = x_k - \frac{f'(x_k)}{f''(x_k)} \] 于是，一直这样迭代下去就可以得到我们的极值了。当然，我们也可以看出来牛顿法的局限性，那就是：如果牛顿法的优化目标是一个复杂的函数，那么我们很难手推出它的一阶导数和二阶导数，甚至有时候，函数可能并不存在一阶导数和二阶导数。

当然了，这并不影响我们写个程序实现它：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import numpy as np
from sympy import *

# 使用牛顿法去找一个一维函数的极小值
def Newton(aimFunc, verb, stPoint = None, step = 10):
    """
    使用牛顿法确定一个函数的最小值
    @param aimFunc 目标函数，sympy对象
    @param verb    变量，也是求导对象
    @param stPoint 初始位置，如果没有就随机一个
    @param step    迭代次数
    """

    if(stPoint is None):
        stPoint = ( np.random.random() - 0.5 ) * 1e5
    curPoint = stPoint
    for i in range(step):
        # 一阶导数值
        onP = diff(aimFunc, verb).evalf(subs ={'x1':curPoint})
        # 二阶导数值
        twP = diff(diff(aimFunc, verb), verb).evalf(subs ={'x1':curPoint})
        curPoint -= onP / twP
    return curPoint


x = symbols('x1')
# 定义目标函数
aimFunction = x ** 2 + 2 * x + 1
print(Newton(aimFunction, x))

我们的目标函数是 \[ f(x) = (x + 1 )^2 \] 即使我们把初始值定在了\(1*10^9\)，但是牛顿法依然能在10次迭代后算到极值-1.00000000000000，说明牛顿发收敛很快。

多维牛顿法

很显然，现实世界中，基本上没有什么需要使用nb优化算法的一维问题，现实世界中的问题一般维度很高，在这一小节，我们介绍二维的牛顿法，以此为例，你可以把它推广到\(n\)维。

首先看，对于一维函数，他有导数来代表因变量随自变量变换的趋势。那么多维函数有什么？没错，就是梯度，那么我们就可以将泰勒展开推广到\(n\)维的函数上去了： \[ \varphi(x) = f(x_k)+\nabla f(x_k)(x - x_k) + \frac{1}{2}(x - x_k)^T \nabla^2f(x_k)(x-x_k) \] 其中，\(\nabla f\)称为\(f\)的梯度向量，\(\nabla^2 f\)称为\(f\)的海森矩阵，他们的定义是： \[ \nabla f= \left[\begin{matrix}\frac{\partial f}{\partial x_1} \\\frac{\partial f}{\partial x_2} \\... \\\frac{\partial f}{\partial x_n}\end{matrix}\right],\nabla^2 f= \left[\begin{matrix}\frac{\partial f}{\partial x_1\partial x_1} & \frac{\partial f}{\partial x_1\partial x_2} & ...&\frac{\partial f}{\partial x_1\partial x_n}\\\frac{\partial f}{\partial x_2\partial x_1} & \frac{\partial f}{\partial x_2\partial x_2} & ...&\frac{\partial f}{\partial x_2\partial x_n}\\... & ... & & ...\\\frac{\partial f}{\partial x_n\partial x_1} & \frac{\partial f}{\partial x_n\partial x_2} & ...&\frac{\partial f}{\partial x_n\partial x_n}\end{matrix}\right] \] 那么同理，需要求极值，所以需要是驻点，所以: \[ \nabla\varphi(x) = 0 \] 解得： \[ x_{k+1} = x_k - {\nabla^2 f}^{-1}*\nabla f \] 最后，迭代就完事了。

下面来看代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import numpy as np
from sympy import *

# 使用牛顿法去找一个一维函数的极小值
def NDimNewton(aimFunc, verbs, stPoint = None, step = 10):
    """
    使用牛顿法确定一个函数的最小值
    @param aimFunc 目标函数，sympy对象
    @param verbs   变量，也是求导对象
    @param stPoint 初始位置，如果没有就随机一个
    @param step    迭代次数
    """

    if(stPoint is None):
        stPoint = ( np.random.random(verbs.shape) - 0.5 ) * 1e5
    curPoint = stPoint
    # 一阶梯度向量
    alpha = np.zeros(curPoint.shape[0])
    # 黑塞矩阵
    H     = np.zeros((curPoint.shape[0], curPoint.shape[0]))

    # 根据迭代次数进行迭代
    for i in range(step):
        print(curPoint)
        # 计算偏导数梯度向量
        for j in range(curPoint.shape[0]):
            alpha[j] = diff(aimFunc, verbs[j]).evalf(subs ={'x1':curPoint[0], 'x2':curPoint[1]})

        # 计算黑塞矩阵
        for j in range(curPoint.shape[0]):
            for k in range(curPoint.shape[0]):
                H[j][k] = diff(diff(aimFunc, verbs[j]), verbs[k]).evalf(subs ={'x1':curPoint[0], 'x2':curPoint[1]})
        
    
        print(H)
        # 更新点
        curPoint -= np.dot(np.array(np.linalg.inv(H), dtype=np.float), np.array(alpha, dtype=np.float))
    return curPoint


x1, x2 = symbols('x1 x2')
# 定义目标函数
aimFunction = 60 - 10 * x1 - 4 * x2 + x1 ** 2 + x2 ** 2 - x1 * x2
print(NDimNewton(aimFunction, np.array([x1, x2]), np.array([1e9, 1e9])))

可以看到，下降路径中，一次就下降到最终结果了。。。那我实验报告怎么写呢？？？