最优化技术——牛顿法

最优化技术——牛顿法

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

参考目录

• 牛顿法与拟牛顿法学习笔记（一）牛顿法

解决的问题

解决的是无约束的极小化问题，用数学语言描述如下：
$$\underset{x}{\min }f(x)\text{\hspace{1em}(0)}$$

简单一维牛顿法

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

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

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维的函数上去了：
$$\phi (x)=f(x_k)+\nabla f(x_k)(x-x_k)+\frac{1}{2}(x-x_k{)}^T{\nabla }^{2}f(x_k)(x-x_k)$$
其中，$\nabla f$称为f的梯度向量，${\nabla }^{2}f$称为f的海森矩阵，他们的定义是：
$$\nabla f=\left[\begin{array}{c}\frac{\partial f}{\partial x_1}\\ frac\partial f\partial x_2\\ \dots \\ frac\partial f\partial x_n\end{array}\right],{\nabla }^{2}f=\left[\begin{array}{cccc}\frac{\partial f}{\partial x_1\partial x_1}& \frac{\partial f}{\partial x_1\partial x_2}& \dots & \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}& \dots & \frac{\partial f}{\partial x_2\partial x_n}\\ \dots & \dots & & \dots \\ frac\partial f\partial x_n\partial x_1& \frac{\partial f}{\partial x_n\partial x_2}& \dots & \frac{\partial f}{\partial x_n\partial x_n}\end{array}\right]$$
那么同理，需要求极值，所以需要是驻点，所以:
$$\nabla \phi (x)=0$$
解得：
$$x_{k+1}=x_k-{{\nabla }^{2}f}^{-1}\ast \nabla f$$
最后，迭代就完事了。

下面来看代码：

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])))

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