人工智能之数学基础篇—高等数学基础（下篇）

  最后，来重点介绍一下方向导数和梯度。因为梯度下降法是机器学习领域非常重要的算法之一，也是应用广泛的优化算法之一。在本篇文章中，我将综合实例来重点介绍梯度下降法及其应用，并利用Python语言编程实现。查看本文之前，可以先阅读人工智能之数学基础篇—高等数学基础上篇和人工智能之数学基础篇—高等数学基础中篇。

6 方向导数

6.1 方向导数的定义

  由二元函数偏导数的概念，我们知道函数 $f(x,y)$ 在点 $\left(x_0,y_0\right)$ 处的两个偏导数，分别是函数 $f$ 在点 $\left(x_0,y_0\right)$ 处沿与 $x$ 轴和 $y$ 轴平行的方向的变化率。但是在许多实际问题中，仅仅研究函数沿这两个特殊方向的变化率是远远不够的，还需要研究函数沿各个不同方向的变化率。例如，在大气气象中，就需要研究温度、气压沿不同方向的变化率。这就是所谓的方向导数。

  从上述的分析中可以看出，偏导数反映的是函数沿坐标轴方向的变化率，方向导数本质上研究的是函数在某点处沿某特定方向的变化率问题。

  比如 $z=f(x,y)$ 在点 $p\left(x_0,y_0\right)$ 沿方向 $l$ 的变化率。假设方向如图10所示，$l$ 为 $xOy$ 平面上以 $p$ 为起始点的一条射线，$p^{{}^{\prime }}$ 为方向 $l$ 上的另一点。

图10 函数沿方向 $l$ 的变化率

  由图10可知，$p$ 与 $p^{{}^{\prime }}$之间的距离 $\left|p{p}^{{}^{\prime }}\right|=\rho =\sqrt{\Delta x^2+\Delta y^2}$。

  函数的增量 $\Delta z=f\left(x_0+\Delta x,y_0+\Delta y\right)-f\left(x_0,y_0\right)$。

  考虑函数的增量与这两点间距离的比值，当 $p^{{}^{\prime }}$ 沿着方向 $l$ 趋于 $p$ 时，如果这个比的极限存在，则称这个极限为函数 $f(x,y)$ 在点 $p$ 沿方向 $l$ 的方向导数，记作 $\frac{\partial f}{\partial l}$，即
$$\frac{\partial f}{\partial l}=\underset{\rho \to 0}{\lim }\frac{f\left(x_0+\Delta x,y_0+\Delta y\right)-f\left(x_0,y_0\right)}{\rho }$$

  从定义可知，当函数 $f(x,y)$ 在点 $p\left(x_0,y_0\right)$ 的偏导数 $f_x\left(x_0,y_0\right)$、$f_y\left(x_0,y_0\right)$ 存在时，函数 $f(x,y)$ 在 $p$ 点沿着 $x$ 轴正向单位向量 e 1 = { 1 , 0 } e_{1}=\left \{ 1,0 \right \} e1​={1,0}，$y$ 轴正向单位向量 e 2 = { 0 , 1 } e_{2}=\left \{ 0,1 \right \} e2​={0,1} 的方向导数存在，且其值依次为$f_x\left(x_0,y_0\right)$，$f_y\left(x_0,y_0\right)$；函数 $f(x,y)$ 在 $p$ 点沿着 $x$ 轴负向单位向量 e 1 ′ = { − 1 , 0 } e_{1}^{'}=\left \{ -1,0 \right \} e1′​={−1,0}，$y$ 轴负向单位向量 e 2 ′ = { 0 , − 1 } e_{2}^{'}=\left \{ 0,-1 \right \} e2′​={0,−1} 的方向导数存在，且其值依次为$-f_x\left(x_0,y_0\right)$，$-f_y\left(x_0,y_0\right)$。

  方向导数 $\frac{\partial f}{\partial l}$ 的存在及计算如下定理。

定理5 如果函数 $z=f(x,y)$ 在点 $p\left(x_0,y_0\right)$ 是可微分的，那么函数在该点沿任意方向 $l$ 的方向导数都存在，且有
$$\frac{\partial f}{\partial l}=\frac{\partial f}{\partial x}cos\phi +\frac{\partial f}{\partial y}sin\phi$$
其中 $\phi$ 为 $x$ 轴到 $l$ 方向的转角。

6.2 方向导数的几何意义

  方向导数的几何意义如图11所示，函数 $z=f(x,y)$ 的变化方向为 $l$，方向导数 ${\begin{array}{c}\frac{\partial f}{\partial l}\end{array}|}_p$ 是函数 $z=f(x,y)$ 在点 $p$ 处沿方向 $l$ 的变化率，从方向 $l$ 作垂直于 $xOy$ 平面的一个平面，与曲面相交成一条曲线 $MQ$，${\begin{array}{c}\frac{\partial f}{\partial l}\end{array}|}_p$ 即为曲线 $MQ$ 在 $M$ 点的切线 $MN$ 的斜率。

图11 方向导数的几何意义

  【例13】求函数 $z=x{e}^{2y}$ 在点 $P\left(1,0\right)$ 到点 $Q\left(2,-1\right)$ 方向的方向导数。

  解：这里方向 $l$ 即向量 P Q → = { 1 , − 1 } \overrightarrow{PQ}=\left \{ 1,-1 \right \} PQ ​={1,−1} 的方向，因此 $x$ 轴到方向 $l$ 的转角 $\phi =-\frac{\pi }{4}$，
   

  由方向导数的意义可知，方向导数 ${\begin{array}{c}\frac{\partial f}{\partial l}\end{array}|}_p$ 是 $z=f(x,y)$ 在点 $p$ 处沿方向 $l$ 的变化率，是曲面上过 $M$ 点的一条曲线的切线的斜率，那么点 $p$ 在某一方向的变化率，就是过曲面 $M$ 点的某一条曲线的切线的斜率。根据曲面的切平面与法线的相关概念可知，曲面上通过点 $M$ 且在点 $M$ 处具有切线的曲线，它们在点 $M$ 处的切线都在同一个平面上，这个平面称为切平面。因此方向导数 ${\begin{array}{c}\frac{\partial f}{\partial l}\end{array}|}_p$ 为过 $M$ 点的切平面上的某条直线的斜率，那么哪一个方向的变化率最大呢?下面继续讨论。

7 梯度

定义16 设函数 $z=f(x,y)$ 在平面区域 $D$ 内具有一阶连续偏导数，则对于每一点 $p\left(x,y\right)ϵD$，都可定出一个向量 $\frac{\partial f}{\partial x}\overrightarrow{i}+\frac{\partial f}{\partial y}\overrightarrow{j}$，这个向量称为函数 $z=f(x,y)$ 在点 $p(x,y)$ 的梯度，记作 $gradf(x,y)$，即
$$gradf(x,y)=\frac{\partial f}{\partial x}\overrightarrow{i}+\frac{\partial f}{\partial y}\overrightarrow{j}$$

  设 $\overrightarrow{e_l}=cos\phi \overrightarrow{i}+sin\phi \overrightarrow{j}$ 是与方向 $l$ 同方向的单位向量（其中 $\phi$ 为 $x$ 轴到方向 $l$ 的转角），则由方向导数的计算公式可知：

  

  由公式推导可以看出，方向导数 $\frac{\partial f}{\partial l}$ 等于梯度与 $\overrightarrow{e_l}$ 的乘积，即梯度的模 $\left|gradf\left(x,y\right)\right|$ 在方向 $l$ 上的投影，当方向 $l$ 与梯度的方向一致时，方向导数 $\frac{\partial f}{\partial l}$ 有最大值。所以，沿梯度方向的方向导数达到最大值，也就是说，梯度方向是函数 $z=f(x,y)$ 在这点增长最快的方向。因此，可以得到如下结论。

  函数在某点的梯度是一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值。

  由梯度的定义可知，梯度的模为
  

  一般来说，函数 $z=f(x,y)$ 在几何上表示一个曲面，这个曲面被平面 $z=c$ （$c$ 是常数）所截得的曲线 $l$ 的方程为：
$$\{\begin{array}{ll}z=f(x,y)& \\ z=c& \end{array}$$
这条曲线 $l$ 在 $xOy$ 平面上的投影是平面曲线 $L^{\ast }$，它在 $xOy$ 平面直角坐标系中的方程为 $f(x,y)=c$。曲线 $L^{\ast }$ 上的任意点，函数值都是 $c$，所以称平面曲线 $L^{\ast }$ 为函数 $z=f(x,y)$ 的等高线。

  由于等高线 $f(x,y)=c$ 上任意点 $(x,y)$ 处的法线的斜率为：
$$-\frac{1}{\frac{dy}{dx}}=-\frac{1}{\left(-\frac{f_x}{f_y}\right)}=\frac{f_y}{f_x}$$

  所以梯度
$$gradf(x,y)=\frac{\partial f}{\partial x}\overrightarrow{i}+\frac{\partial f}{\partial y}\overrightarrow{j}$$为等高线上点 $(x,y)$ 处的法向量，因此可得到梯度与等高线的关系如下。

  如图12所示，图中显示了6条等高线，在等高线 $f(x,y)=c$ 上有一点 $p_0\left(x_0,y_0\right)$，可以看出点 $p_0\left(x_0,y_0\right)$ 的切线与法线垂直。图中 $c_2>c>c_1$，$p_0\left(x_0,y_0\right)$ 点的法线方向与梯度 $▽f\left(x_0,y_0\right)$ 相同。

图12 等高线、切线、法线的关系

  梯度为方向导数取最大值的方向的概念并不好理解。再通俗一点，设想曲面为一座高山，而这座山像一个球的半面，如图13所示。现在，有一个人站在半山腰 $p$ 点位置上，他想要寻找下山最快的方向，在他所处的位置做一个切面，所有的方向导数都在此切面上，哪一个方向使他下山最快呢？采用他所处位置的等高线的切线的法向量最快，由于是下山，因此采用的是梯度的反方向。

图13 梯度的通俗解释

  【例14】求函数 $u=xyz+z^2+5$，求 $gradu$ 及在点 $M(0,1,-1)$ 处方向导数的最大值。

  解：

8 综合实例一梯度下降法求函数的最小值

  梯度下降法，是寻找函数极小值最常用的优化方法。当目标函数是凸函数时，梯度下降法的解是全局解。但在一般情况下，其解不保证是全局最优解。最普遍的做法是，在已知参数当前值时，按当前点对应的梯度向量的反方向及事先定好的步长大小对参数进行调整。按上述方法对参数做出多次调整之后，函数就会通近一个极小值。

  例如，假设函数 $f(x)$ 为一元连续函数，初始值为 $x_0$，$\alpha$ 为步长，已知在点 $x_0$ 的梯度为 $gradf\left(x_0\right)$，那么下一个点的坐标为 $x_1=x_0-\alpha \ast gradf\left(x_0\right)$，然后再求点 $x_1$ 的梯度，反复迭代，直到 $f\left(x_n\right)-f\left(x_{n-1}\right)$ 的绝对差极小，迭代结束，此时的 $f\left(x_n\right)$ 即为极小值。

  下面举一个例子来说明梯度下降法的基本原理。

  假设曲面上一只蚂蚁突遇火灾，该如何快速逃跑呢？此问题可以类比为一个下山的过程：假设蚂蚁被困在高山上，需要快速找到山的最低点，即山谷。但此时山上的烟雾很大，可视度很低，下山的路径无法确定，它必须利用自己周围的信息去寻找。这个时候，它就可以利用梯度下降法来帮助自己下山。以它当前所处的位置为基准，寻找这个位置最陡峭的地方，朝着山的高度下降最快的地方走。但是，由于山的地形复杂，坡度变化随机，这样一直往下走无法确定路径是否正确，因此可以采用每走一段距离，都反复采用同一个方法（即寻找当前位置最陡峭的地方，朝着山的高度下降最快的地方走），最后就能成功地抵达山谷。

  梯度下降的过程就与蚂蚁下山的场景类似。

  首先，有一个可微分的函数 $f(x)$，这个函数就代表着一座山，我们的目标就是找到这个函数的最小值，也就是山底。根据之前的场景假设，最快下山的方式就是找到当前位置最陡峭的方向，然后沿着此方向向下走，对应到函数中，就是找到给定点 $x_0ϵI$ 的梯度，梯度相反的方向就能让函数值下降最快。重复利用这个方法，反复求取梯度，最后就能到达局部的最小值。梯度下降法搜索迭代过程如图14所示。

图14 梯度下降法的迭代过程

  【例15】应用 $Python$ 编程实现梯度下降法求解下面函数的最小值，并使用 $Matplotlib$、$mp{l}_{-}toolkits$ 库画出函数的图形。
$$minf\left(x\right)=x_1-x_2+2x_1^2+2x_1{x}_2+x_2^2$$给定初值 $X^{\left(0\right)}={\left(0,0\right)}^T$。

  解：$matplotlib.pyplot$ 是一个有命令风格的函数集合，它看起来和 $MATLAB$ 很相似，每一个 $pyplot$ 函数都使一幅图像 $figure$ 做出些许改变；$mp{l}_{-}toolkits.mplot3d$ 是三维工具包。本例中 $mplot3d$ 仍使用 $figure$ 对象，只不过 $Axes$ 对象要替换为该工具集的 $Axes3d$ 对象。

  【代码如下】

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D


# f(w1,w2) = w1 - w2 + 2*w1^2 + w2^2 + 2*w1*w2
def target_function(w):      # 目标函数
    w_1, w_2 = w
    return w_1 - w_2 + 2 * w_1 ** 2 + w_2 ** 2 + 2 * w_1 * w_2


def gradient_function(w):    # 梯度函数：分别对w1,w2求偏导
    w_1, w_2 = w
    w1_grad = 1 + 4 * w_1 + 2 * w_2
    w2_grad = -1 + 2 * w_1 + 2 * w_2
    return np.array([w1_grad, w2_grad])


# 梯度下降算法
def batch_gradient_distance(target_func, gradient_func, init_w, learning_rate=0.0008, tolerance=7e-9):
    w = init_w
    target_value = target_func(w)
    counts = 0  # 用于计算次数
    while counts < 50000:
        gradient = gradient_func(w)
        next_w = w - gradient * learning_rate
        next_target_value = target_func(next_w)
        if abs(next_target_value-target_value) < tolerance:
            print("此结果经过了", counts, "次循环")
            return next_w
        else:
            w, target_value = next_w, next_target_value
            counts += 1
    else:
        print("没有取到极值点")


if __name__ == '__main__':
    init_w = np.array([0, 0])   # 初始坐标点为（0，0）
    w1, w2 = batch_gradient_distance(target_function, gradient_function, init_w)
    print(w1, w2)
    # 画图
    x1 = np.arange(-2, 2, 0.1)
    x2 = np.arange(-2, 2, 0.1)

    x1, x2 = np.meshgrid(x1, x2)    # meshgrid :3D坐标系

    z = x1 - x2 + 2 * x1 ** 2 + x2 ** 2 + 2 * x1 * x2

    fig = plt.figure()              # figure 对象
    ax = Axes3D(fig, auto_add_to_figure=False)  # Axes3D 对象
    fig.add_axes(ax)
    # fig = plt.figure()
    # ax = fig.gca(projection='3d')

    ax.plot_surface(x1, x2, z, rstride=1, cstride=1, cmap="rainbow")    # 绘制3D坐标系中的函数图像
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')

    ax.scatter(w1, w2, target_function([w1, w2]), s=30, c='red')        # 绘制已经找到的极值点
    ax.text(-2, -1, 10, '(' + str(w1) + "," + str(w2) + ')', color='red')
    plt.show()

  【运行结果】

图15 函数的图形

9 Python中相关库简介

  在高等数学基础篇的三篇文章中利用了Python来求解相关的问题，主要用到了SymPy 、 NumPy 、 SciPy 、 Matplotlib 库 。因此，接下来将对几个相关的库做简要总结。

9.1 SymPy库简介

  SymPy 是 Python 的一个科学计算库，用强大的符号计算体系完成诸如多项式求值、求极限 、求导 、解方程 、求积分 、解微分方程、级数展开 、矩阵运算等计算 。

  （1）常用的 SymPy 内置符号

  自然对数的底 $e$ 的表示方式 。

  【代码如下】

import sympy

print(sympy.E)
print(sympy.log(sympy.E))       # 自然对数的底 e 的表示方式。

E
1

  无穷大 $\infty$ 的表示方式 。

  【代码如下】

import sympy

print(1/sympy.oo)               # 无穷大 ∞ 的表示方式。

0

  圆周率 $\pi$ 的表示方式 。

  【代码如下】

import sympy

print(sympy.sin(sympy.pi))      # 圆周率 π 的表示方式。

0

  （2）用 SymPy 进行初等运算

  常用的函数有：求对数函数 sympy.log 、正弦函数 sympy.sin 、求平方根函数 sympy.sqrt 、求 n 次方根函数 sympy.root、求阶乘函数 sympy.factorial 等。

  （3）表达式与表达式求值

  SymPy 可以用一套符号系统来表示一个表达式，如函数、多项式等，并且可以进行求值。

  【代码如下】

import sympy

x = sympy.Symbol('x')
fx = 2*x + 1    # fx 是一个表达式
print(fx.evalf(subs={x: 2}))        # 用 evalf() 函数，传入变量的值，对表达式进行求值

x, y = sympy.symbols('x y')
f = 2 * x + y
f.evalf(subs={x: 1, y: 2})          # 以字典的形式传入多个变量的值

print(f.evalf(subs={x: 1}))         # 如果只传入一个变量的值，则输出原来的表达式
print(f.evalf(subs={x: 1, y: 2}))   # 用 evalf() 函数，传入变量的值，对表达式进行求值

5.00000000000000
y + 2.0
4.00000000000000

  （4）利用 SymPy 求极限

  使用函数 sympy.limit 求极限，可以参考人工智能之数学基础篇—高等数学基础（上篇）中的例5、例6和例8。

  （5）利用 SymPy 求导

  使用函数 sympy.diff 求导数，传入两个参数即函数表达式和变量名，可以参考人工智能之数学基础篇—高等数学基础（中篇）中的例10和例11。

9.2 NumPy库简介

  NumPy 是 Python 的一个扩展程序库，支将大维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。NumPy 是一个运行速度非常快的数学库，主要用于数组计算。

  NumPy 通常与 SciPy 和 Matplotlib （绘图库）一起使用， 这种组合广泛应用于替代MATLAB，是一个强大的科学计算环境。 在本篇文章例15中就利用了这些库实现了求函数最小值的解，并用 NumPy 和 Matplotlib 的一些函数实现了函数图像的可视化。

  （1）数组的操作

  【代码如下】

import numpy as np

a = np.array([1, 2, 3])                 # 创建一个秩 (rank)1 数组
print(a[0], a[1], a[2])
b = np.array([[1, 2, 3], [4, 5, 6]])    # 创建一个秩 2 数组
print(b[0, 0], b[0, 1], b[1, 0])

1 2 3
1 2 4

  （2）在绘制三维图表时，需要用到 NumPy 中的 mgrid 函数。它会返回一个密集的多维网格，一般形式为 np.mgrid[start: end: step]，其中 start 表示开始值，end 表示结束坐标（不包含此点）， step表示步数。

  【代码如下】

import numpy as np

print(np.mgrid[-1:4:2])

[-1  1  3]

print(np.mgrid[-1:4:2, -3:1:1])

[[[-1 -1 -1 -1]
  [ 1  1  1  1]
  [ 3  3  3  3]]

 [[-3 -2 -1  0]
  [-3 -2 -1  0]
  [-3 -2 -1  0]]]
 
print(np.mgrid[-1:4:2, -3:1:2])

[[[-1 -1]
  [ 1  1]
  [ 3  3]]

 [[-3 -1]
  [-3 -1]
  [-3 -1]]]

9.3 SciPy库简介

  SciPy 是世界上著名的 Python 开源科学计算库，构建在 NumPy 之上，功能更为强大。SciPy 函数库在 NumPy 库的基础上增加了众多的数以及工程计算中常用的库函数。可以说，NumPy 是一个纯数学的计算模块， 而 SciPy 是一个更高阶的科学计算库 。SciPy 库需要 NumPy 库的支持。由于这种依赖关系，NumPy 库的安装要先于 SciPy 库。后面会使用该库来计算函数的积分。

9.4 求导的3种方式

  （1）使用 SymPy 的 diff 函数，可以得到 $f(x)$ 的导数表达式；

  （2）使用 scipy.misc 模块下的 derivative 函数。SciPy 的求导相对简单也容易理解。

  【代码如下】

import numpy as np
from scipy.misc import derivative

def f(x):
    return x**5
print(derivative(f, 2, dx=1e-6))   # 对函数在 x=2 处求导

80.00000000230045

  （3）使用 NumPy 模块里 polyld 函数构造 $f(x)$，polyld 函数的形参是多项式的系数，最左侧的是最高次数的系数，构造的函数为多项式。NumPy 的 polyder 函数和 deriv 函数的作用差不多，都是对多项式求导，可以得到函数导数的表达式和在某点的导数。

  【例16】对多项式 $x^5+2x^4+3x^2+5$ 求导。

  【代码如下】

import numpy as np

p = np.poly1d([1, 2, 0, 3, 0, 5])   # 构造多项式
print(p)  # 下面两行为多项式的显示形式，5、4、2 是下一行项数所对应的幂次。

   5     4     2
1 x + 2 x + 3 x + 5

print(np.polyder(p, 1))       # 求一阶导数

   4     3
5 x + 8 x + 6 x

print(np.polyder(p, 1)(1.0))  # 求一阶导数在点 x=1 的值。

19.0

print(p.deriv(1))             # 求一阶导数

   4     3
5 x + 8 x + 6 x

print(p.deriv(1)(1.0))        # 求一阶导数在点 x=1 的值。

19.0