数值计算之 梯度下降法与函数极值

前言

本篇将开始介绍优化算法。首先是梯度下降法，在最小二乘与深度学习中，都是最常用的最优化求解方法和思想。

微积分基础

一元函数的极值，导数与泰勒展开

对于一元函数$f(x)$而言，当$x_0$满足以下条件时，$f(x_0)$取得极值：
$$f^{\prime }(x_0)=0,f^{\prime \prime }(x_0)\ne 0$$
将$f(x)$在$x_0$处泰勒展开：
$$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2!}{f}^{\prime \prime }(x_0)(x-x_0{)}^2+o((x-x_0{)}^2)$$
可以从泰勒展开中看出，当$f^{\prime }(x)=0,f^{\prime \prime }(x_0)>0$，在$x_0$附近必然有$f(x)>f(x_0)$成立；当$f^{\prime }(x)=0,f^{\prime \prime }(x_0)<0$，在$x_0$附近必然有$f(x)<f(x_0)$成立。这就把一元函数的局部最值与泰勒公式联系起来了。

多元函数的泰勒展开

对于实值向量函数$f(\mathbf{x})$，其在${\mathbf{x}}_0$处的泰勒展开可表示为：
$$f(\mathbf{x})=f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0)+(\mathbf{x}-{\mathbf{x}}_0{)}^{T}H({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)+o^n$$
类似的，当梯度$\nabla f({\mathbf{x}}_0)=0$, 海森矩阵$H$正定时，$f({\mathbf{x}}_0)$是一个局部极值点。

从多元函数的泰勒展开也能看出，如果${\mathbf{x}}_0$处的梯度不为$0$，则在周围存在某一点$\mathbf{x}$使得$f(\mathbf{x})>f({\mathbf{x}}_0)$，同时也存在某一点$\mathbf{x}$使得$f(\mathbf{x})<f({\mathbf{x}}_0)$。

自然而然的想到，如果我们要寻找函数的极大值，可以从${\mathbf{x}}_0$开始，选择周围的某一点$\mathbf{x}=\Delta \mathbf{x}+{\mathbf{x}}_0$，使得$f(\mathbf{x})>f({\mathbf{x}}_0)$，到达点${\mathbf{x}}_1$，然后再选择周围某一个点$\mathbf{x}=\Delta \mathbf{x}+{\mathbf{x}}_1$，继续迭代到满足局部极值条件为止。寻找极小值同理。

这就是梯度下降（上升）法的思想。

梯度下降法

以上迭代具有两个核心问题：①如何选择周围点，也就是如何选择$\Delta \mathbf{x}$；②如何判断满足局部极值条件，也是就是什么时候结束迭代。下面以最常见的寻找局部最小为例。

梯度方向

首先讨论问题①，回到多元泰勒展开式：
$$f(\mathbf{x})=f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_0{)}^{T}H({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)+o^n$$
梯度$\nabla f({\mathbf{x}}_0)\ne 0$。假设我们要寻找的$\Delta \mathbf{x}$的长度固定，则每次迭代的函数增量为：
$$\Delta f=f(\mathbf{x})-f({\mathbf{x}}_0)=\nabla f({\mathbf{x}}_0)\cdot (\Delta \mathbf{x})+{\mathbf{o}}^{\mathbf{n}}$$
等式右边是一个内积，可以简化表示为：
$$\begin{gathered} g_{{\mathbf{x}}_0}=\nabla f({\mathbf{x}}_0) \\ {e}_{\Delta \mathbf{x}}=\Delta \mathbf{x} \\ \Delta f=f(\mathbf{x})-f({\mathbf{x}}_0)=g_{{\mathbf{x}}_0}\cdot e_{\Delta \mathbf{x}}=|g_{{\mathbf{x}}_0}||e_{\Delta \mathbf{x}}|\cos \theta \end{gathered}$$
其中，$\cos \theta$是梯度向量$\nabla f({\mathbf{x}}_0)$与$\Delta \mathbf{x}$的夹角。$\theta =0$，$\Delta f>0$并且取最大值；$\theta =\pi$，$\Delta f<0$且取最小值。

由此得出一个结论：梯度方向是函数上升最快的方向，梯度反方向是函数下降最快的方向。

因此，在迭代过程中，我们选择的$\Delta \mathbf{x}$的方向与梯度方向相反，即：
$$\Delta \mathbf{x}=-\frac{||\Delta \mathbf{x}|{|}_2}{||\nabla f({\mathbf{x}}_0)|{|}_2}\nabla f({\mathbf{x}}_0)$$

另外，可以通过学习率$\lambda$控制$\Delta \mathbf{x}$的长度。一方面是由于泰勒展开只在${\mathbf{x}}_0$的局部范围才成立，另一个方面是较大的$\Delta \mathbf{x}$可能会导致函数找不到局部极值区域，而较小的$\Delta \mathbf{x}$需要更多的迭代次数。梯度下降法的迭代方程可表示为：
$$h(\mathbf{x})=\mathbf{x}-\lambda \nabla f(\mathbf{x})$$

终止条件

对于问题②，有比较多的方式来终止迭代，比如设置最大迭代次数，设置增量阈值，梯度限制等。

可以设置函数值增量$\Delta f$阈值来控制迭代过程：
$$\begin{gathered} if\Delta f_{last}<threshold \\ then\mathbf{x}=h_{last}(\mathbf{x}) \end{gathered}$$
也可以设置$\Delta \mathbf{x}$的长度阈值：
$$\begin{gathered} if\Delta \mathbf{x}<threshold \\ then\mathbf{x}={\mathbf{x}}_{last} \end{gathered}$$

代码举例

这里给出一个我自己写的求二元函数$f(x,y)=x^2+y^2$与$f(x,y)=x^2+2y^2+2xy+3x-y-2$的极小值的代码：

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 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 grad_decent(x, y, F, lr):
    grad_x, grad_y = partial_derivate_xy(x, y, F)
    x_opt = x - lr * grad_x
    y_opt = y - lr * grad_y
    return x_opt, y_opt


def optimizer(x0, y0, F, lr=0.01, th=0.0001):
    x = x0
    y = y0
    counter = 0
    while True:
        x_opt, y_opt = grad_decent(x, y, F, lr)
        if (x_opt - x)**2 + (y_opt - y)**2 < th**2:
            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)
    print(result_x, result_y)

后记

本篇记录了梯度下降法的原理，后续会再总结SGD，Adam等优化器。