数值计算之 牛顿法与函数极值

前言

本篇继续优化理论的算法学习，牛顿法。

最速下降法

首先回顾上次提到的梯度下降法（其实就是最速下降法）：通过求取多元函数在某个点处的梯度，沿着梯度的反方向前进，直到达到迭代停止条件。

对于多元函数（实值向量函数）$f(\mathbf{x})$，其在${\mathbf{x}}_0$处的泰勒展开可表示为：
$$f(\mathbf{x})=f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0{)}^T(\mathbf{x}-{\mathbf{x}}_0)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_0{)}^T\mathbf{H}({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)$$
其中，$\nabla ,H$分别是函数梯度，海森矩阵。

梯度下降法直接取一阶梯度的反方向作为优化方向，因此称为最速下降法（每次迭代的方向都是下降最快的方向）。

牛顿法

回到多元泰勒展开：
$$f(\mathbf{x})=f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_0{)}^T\mathbf{H}({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)$$
$f(\mathbf{x})$的极小值处的导数应当等于$0$。现在取初始点${\mathbf{x}}_0$，将$f(\mathbf{x})$在${\mathbf{x}}_0$处的展开式对$\mathbf{x}$进行求导：
$$\begin{gathered} \nabla f({\mathbf{x}}_0{)}^T+{(\mathbf{x}-{\mathbf{x}}_0)}^T\mathbf{H}({\mathbf{x}}_0)=0 \\ \nabla f({\mathbf{x}}_0)=-\mathbf{H}({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0) \\ \mathbf{x}={\mathbf{x}}_0-{\mathbf{H}({\mathbf{x}}_0)}^{-1}\nabla f({\mathbf{x}}_0) \end{gathered}$$
这样就得到了新的$\mathbf{x}$，这就是牛顿法。

可以给出更具体的牛顿法求极值过程：

1. 确定实值向量函数$f(\mathbf{x})$的初始点${\mathbf{x}}_0$，梯度$\nabla f$的表达式，海森矩阵$\mathbf{H}(\mathbf{x})$的表达式，终止条件$\delta ,ϵ$
2. 计算梯度$\nabla f({\mathbf{x}}_0)$，如果$||\nabla f({\mathbf{x}}_0)||<\delta$，终止迭代，得到解$\mathbf{x}={\mathbf{x}}_0$；否则进入第3步
3. 计算海森矩阵，并迭代$${\mathbf{x}}^{\prime }={\mathbf{x}}_0-{\mathbf{H}({\mathbf{x}}_0)}^{-1}\nabla f({\mathbf{x}}_0)$$
4. 计算$\Delta f=f({\mathbf{x}}^{\prime })-f({\mathbf{x}}_0)$，如果$\Delta f<ϵ$，终止迭代，得到解$\mathbf{x}={\mathbf{x}}^{\prime }$；否则${\mathbf{x}}_0={\mathbf{x}}^{\prime }$，回到第2步

牛顿法分析

现在考虑牛顿法的优化迭代表达式：
$${\mathbf{x}}^{\prime }={\mathbf{x}}_0-{\mathbf{H}({\mathbf{x}}_0)}^{-1}\nabla f({\mathbf{x}}_0)$$

将牛顿法的迭代方向与梯度方向做内积：
$$\begin{gathered} -{\mathbf{H}({\mathbf{x}}_0)}^{-1}\nabla f({\mathbf{x}}_0)\cdot \nabla f({\mathbf{x}}_0) \\ =-{\nabla f({\mathbf{x}}_0)}^T{\mathbf{H}({\mathbf{x}}_0)}^{-T}\nabla f({\mathbf{x}}_0) \\ =-{\nabla f({\mathbf{x}}_0)}^T{\mathbf{H}({\mathbf{x}}_0)}^{-1}\nabla f({\mathbf{x}}_0) \end{gathered}$$
当海森矩阵是正定矩阵时，迭代方向与梯度方向的内积为负，即函数值在不断减小。这就涉及到了优化理论的重要问题：什么条件能够满足海森矩阵的正定性？

牛顿法的优点：优化收敛速度比梯度下降法更快。原因是牛顿法不仅考虑到函数的一阶梯度，也考虑到了海森矩阵（二阶梯度），在每个优化点处附近都能一次求得最佳迭代值。

牛顿法的缺点：每次迭代的速度较慢。这是因为迭代值的计算中涉及到了海森矩阵的逆，使得单步迭代效率比较低。

代码示例

和上次一样，我写了一个同牛顿法寻找二元函数极小值的代码，二元五次函数，如下所示：

import numpy as np


def partial_derivate_xy(x, y, F):
    dx = (F(x + 0.001, y) - F(x, y))/0.001
    dy = (F(x, y + 0.001) - F(x, y))/0.001
    return dx, dy


def partial_partial_derivate_xy(x, y, F):
    dx, dy = partial_derivate_xy(x, y, F)
    dxx = (partial_derivate_xy(x + 0.001, y, F)[0] - dx) / 0.001
    dyy = (partial_derivate_xy(x, y + 0.001, F)[1] - dy) / 0.001
    dxy = (partial_derivate_xy(x, y + 0.001, F)[0] - dx) / 0.001
    dyx = (partial_derivate_xy(x + 0.001, y, F)[1] - dy) / 0.001
    return dxx, dyy, dxy, dyx


def non_linear_func(x, y):
    fxy = 0.5 * (x ** 2 + y ** 2)
    return fxy


def non_linear_func_2(x, y):
    fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2
    return fxy


def non_linear_func_3(x, y):
    fxy = 0.5 * (x ** 2 - y ** 2)
    return fxy


def non_linear_func_4(x, y):
    fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5
    return fxy


def newton_optim(x, y, F):
    dx, dy = partial_derivate_xy(x, y, F)
    dxx, dyy, dxy, dyx = partial_partial_derivate_xy(x, y, F)
    grad = np.array([[dx], [dy]])
    hessian = np.array([[dxx, dxy], [dyx, dyy]])
    hessian_inv = np.linalg.inv(hessian)
    vec_delta = np.matmul(hessian_inv, grad)
    vec_opt = np.array([[x], [y]]) - vec_delta
    x_opt = vec_opt[0][0]
    y_opt = vec_opt[1][0]
    return x_opt, y_opt, vec_delta


def optimizer(x0, y0, F, th=0.01):
    x = x0
    y = y0
    counter = 0
    while True:
        x_opt, y_opt, vec_delta = newton_optim(x, y, F)
        if np.linalg.norm(vec_delta) < th:
            break
        x = x_opt
        y = y_opt
        counter = counter + 1
        print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y))
    return x, y


if __name__ == '__main__':
    x0 = 2.
    y0 = 2.
    result_x, result_y = optimizer(x0, y0, non_linear_func_4)
    print(result_x, result_y)

牛顿法迭代11次后得到优化解。

梯度下降法在lr=0.01时迭代149次得到优化解。

运行1000次时，牛顿法用时0.586s，梯度下降法用时0.285s，可见虽然牛顿法迭代次数少，但速度慢。

后记

下篇将进入高斯牛顿法，后续还会有LM算法、拟牛顿法等。