数值求导和自动求导

本文转载自查看原文 2021-03-04 15:48 398 人工智能

2021-03-04

数值求导和自动求导

早在高中阶段，我们就开始接触导数，了解过常用函数的求导公式。大学时，我们进一步懂得了用极限定义导数，比如，函数 $f$ 在 $x$ 处的导数定义为

$f'\left( x \right) = \lim_{\epsilon \rightarrow 0}{\frac{f\left( x + \epsilon \right) - f \left( x \right)}{\epsilon}} \\$

然而，这个定义式似乎从来没有派上过用场，始终束之高阁。因为对我们来说，这个式子是没法计算的， $\epsilon$ 趋近于 $0$ ，超出了我们手工计算的能力范畴。另一方面，各种求导公式和求导法则可以让我们直接求出导函数的解析解，所以完全不需要根据定义来计算某个点的导数。

但出了校园，形势就变得微妙起来。很多时候，我们无法求得导函数的解析解，或者求解析解的代价太大。这种情况下，数值求导和自动求导应运而生，这是数学与计算机科学碰撞出的火花，我们将会看到它在工程领域绽放的光芒。

数值求导

如果函数 $f$ 是个黑盒子，即我们不知道它的解析形式，但给定任意的输入 $x$ ，我们都能得到唯一的输出 $f(x)$ 。这就是典型的计算机程序的特点，函数的内部实现被封装在库中，调用方得不到库的源码，这时候要想求函数的导数，就只能用数值求导方法。

本文介绍的数值求导方法，也是最常用的数值求导方法，称为有限差分（Finite-Difference）。该方法正是源自文章开头介绍的导数定义式。应用该式时，我们需要去掉 $\lim_{\epsilon \rightarrow 0}{}$ 运算符，用一个实际的小量 $\epsilon$ 代入，但是这样做势必会引入误差。显然， $\epsilon$ 越小，误差也会越小，但到底误差与 $\epsilon$ 之间是怎样的关系，线性关系还是二次关系，就需要严谨的推导才能得知。

使用泰勒定理展开 $f(x+\epsilon)$ 可得

$f(x+\epsilon) = f(x) + f'(x)\epsilon + \frac{1}{2}f''(x + t\epsilon)\epsilon^2, \quad \text{some } t \in (0, 1) \tag{7.2} \\$

如果函数的二阶导有限，即存在使得 $\| f''(x) \| \le L$ 成立的 $L$ ，我们就可以把式(7.2)做一次放缩，得到

$\| f(x+\epsilon) - f(x) - f'(x)\epsilon \| \le (L/2) \epsilon^2 \tag{7.3} \\$

进而求得导数

$f'(x) = \frac{f(x + \epsilon) - f(x)}{\epsilon} + \delta_\epsilon, \quad \text{where } | \delta_\epsilon| \le (L/2)\epsilon \tag{7.4} \\$

式(7.4)把式(7.3)中的不等号转化为一个小量 $\delta_\epsilon$ ，该小量的上界由 $\epsilon$ 决定。此时，我们可以说，有限差分的结果与真实导数值之间的误差在 $O(\epsilon)$ 这个量级。

看起来似乎让 $\epsilon$ 越小越好，但实际并没有这么简单。计算机的精度是有上限的，任何运算都不可能得到完全精确的结果，而是会存在舍入误差。假设用来计算式 $(7.4)$ 的量都用双精度浮点数表示，我们来看看计算机中这些数会有怎样的舍入误差。在学习C语言的时候一般会提到，双精度浮点数由64个二进制位组成，IEEE规定它们的排布如下

最高位是符号位，接下来的11位是指数位，剩下的52位是小数位。现在我们考虑一个问题，对于任意一个可以用浮点数表示的实数，它的真实值与其浮点数表示的值之间相差多少？咋看上去，这个问题的答案似乎是浮点数所能表示的最接近0的实数，也就是指数位取0、小数位取1对应的数。但实际上，由于有效数字位数只有52位，指数越大，有效数字的最低位对应的大小就越大，所以答案是不确定的。我们只能确保真实值与其浮点数表示的值的前52位二进制位是完全相同的，换算成10进制大概是15位。举例来说， $10^{15}$ 与其浮点数表示相差大约1，而1与其浮点数表示相差大约 $10^{-15}$ 。由此推广到任意一个数 $x$ ，它的浮点数可以表示为

$\text{fl}(x)=x(1+\epsilon), \quad \text { where }|\epsilon| \leq \mathbf{u} \tag{A.58} \\$

其中 $\mathbf{u} \approx 10^{-15}$ 。可以发现， $\mathbf{u}$ 其实表示的是相对精度。现在回到式(7.4)，如果我们用计算机计算的结果代替理想值，即用 $\text{fl}(f(x+\epsilon) )$ 替代 $f(x+\epsilon)$ ，用 $\text{fl}(f(x))$ 替代 $f(x)$ ，结果势必会变化，由式(A.58)可知

$\begin{aligned} | \text{fl}(f(x)) - f(x)| &= |f(x) + f(x)\epsilon - f(x)| \\ &= |f(x) \epsilon| \\ &\le L_f |\epsilon| \\ &\le L_f \mathbf{u} \end{aligned} \\$

其中， $L_f$ 是函数值的界限。同理，

$| \text{fl}(f(x+\epsilon)) - f(x+\epsilon)| \le L_f \mathbf{u} \\$

这意味着

$\text{fl}(f(x)) = f(x) + \delta_{\epsilon_1}, \quad \text{where } |\delta_{\epsilon_1}| \le L_f \mathbf{u} \\ \text{fl}(f(x + \epsilon)) = f(x + \epsilon) + \delta_{\epsilon_2}, \quad \text{where } |\delta_{\epsilon_2}| \le L_f \mathbf{u} \\$

此时，我们再按照式(7.4)的计算方式计算导数的浮点数表示

$\begin{aligned} \text{fl}(f'(x)) &= \frac{\text{fl}(f(x + \epsilon)) - \text{fl}(f(x))}{\epsilon} + \delta_\epsilon\\ &= \frac{f(x + \epsilon) - f(x) + \delta_{\epsilon_2} - \delta_{\epsilon_1}}{\epsilon} + \delta_\epsilon \\ &= \frac{f(x + \epsilon) - f(x)}{\epsilon} + \delta \end{aligned} , \quad \text{where } | \delta| \le (L/2)\epsilon + 2L_f\mathbf{u}/\epsilon \tag{7.5} \\$

现在，误差不只和 $\epsilon$ 有关，还与计算机的相对精度 $\mathbf{u}$ 有关。如果我们把 $\delta$ 对 $\epsilon$ 求导，令导数等于0，可以求得当 $\epsilon^2 = \frac{4L_f\mathbf{u}}{L}$ 时， $\delta$ 取极小值。这个结论说明， $\epsilon$ 并不是越小越好。实际应用中，我们假设 $\frac{4L_f}{L}$ 不会太大也不会太小，所以取 $\epsilon = \sqrt{\mathbf{u}}$ 即可。

这种有限差分方法称为前向差分（forward-difference），因为只计算了当前点和前面的点的函数值，可以想象，这种方法算出来的结果总是会有一些偏差。更好的方法是中心差分（central-difference），定义如下

$f'(x) \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2\epsilon} \tag{7.7} \\$

如果我们仍按照前面的方式分别对 $f(x+\epsilon)$ 和 $f(x-\epsilon)$ 泰勒展开，可以发现式(7.7)的近似误差在 $O(\epsilon^2)$ 数量级。

那是不是中心差分一定比前向差分好呢？未必。实际中，函数的自变量 $x$ 通常是一个向量。此时导数也是个向量，有限差分方法需要分别对自变量的各个分量求导。假设 $x$ 是 $N$ 维向量，那么前向差分需要计算 $N+1$ 次函数值，而中心差分则需要计算 $2N$ 次函数值。可见，中心差分精度高的代价是计算量大，实际使用时需要有所取舍。

稀疏雅克比

目标函数的导数又称为雅克比矩阵。对于标量函数来说，雅可比矩阵是个行向量，形式如下

$\nabla f(x)=\left[\begin{array}{} \frac{\partial f}{\partial x_1} & \frac{\partial f}{\partial x_2} & \cdots & \frac{\partial f}{\partial x_n} \end{array}\right] \\$

刚刚提到过，这种函数的前向差分需要计算 $N+1$ 次函数值。而对于向量函数，雅可比矩阵是个 $M×N$ 的矩阵， $M$ 是目标函数的维度，形式如下

$\begin{equation*} \nabla f(x) = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots &\frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n} \end{bmatrix} \end{equation*} \\$

这时候，矩阵的每一个元素都需要计算一次函数值（这里所说的函数值是指向量函数中的结果向量的每一维），总共需要计算 $M(N+1)$ 次。

于是，输入变量和目标函数的维度越高，雅可比的计算量越大，数值求导很容易变得不再可行。好在某些实际问题中，雅可比矩阵会呈现出稀疏的特征，比如下面这个样子

$\begin{bmatrix} × & × & & & & \\ × & × & ×& & & \\ & ×& × & × & & \\ & & ×& ×& ×& \\ & & & ×& ×& × \\ & & & & × & × \end{bmatrix} \tag{7.13} \\$

只有对角的条带状元素有值，其它位置都是0。如何利用这一特点加速雅可比计算呢？

把有限微分公式推广到向量情况，对于雅可比矩阵中的任意一个值 $\frac{\partial f_j}{\partial x_i}$ ，需要取 $\epsilon_i = \epsilon e_i$ ，其中 $e_i$ 是第 $i$ 维的单位基向量。仔细观察式(7.13)的第一列和第四列，可以发现， $x_1$ 只对导数中的 $f'_1$ 和 $f'_2$ 分量产生影响， $x_4$ 只对导数中的 $f'_3$ 、 $f'_4$ 和 $f'_5$ 分量产生影响，这说明 $x_1$ 和 $x_4$ 关于导数值是相互独立的。这样我们就可以取 $\epsilon_{1,4} = \epsilon ( e_1 + e_4)$ ，一次性计算出 $\frac{\partial f_j}{\partial x_1}$ 和 $\frac{\partial f_j}{\partial x_2}$ 。

这就是稀疏雅可比带来的效率提升。推广到一般情况，即使雅可比矩阵的形式不是条带状，也存在一些方法可以找到相互独立的若干个分量，从而加速计算。不过可以预见的是，对任意形状的雅可比矩阵搜索独立分量并不容易，它本质上是一个图指派问题（graph assignment）。

条带状的雅可比矩阵并不罕见，如果你接触过视觉SLAM，就会发现视觉SLAM中的雅可比矩阵都是条带状的。这是因为在整个时间序列中，某个时间点的误差项只和少数的若干个位姿和路标点发生关联，你不可能在任何位置观测到所有的路标点。

既然雅可比矩阵可以有限差分，海森矩阵也可以。而且海森矩阵是对称的，所以我们只需要计算一半的元素即可。海森矩阵也有常见的稀疏形式，通常是箭头状，利用这一稀疏性也可以大大减少计算量。由于篇幅限制，本文就不详细介绍了。

自动求导

数值求导很方便，不需要知道目标函数的表达式就可以计算。但是大多数情况是，我们知道目标函数的表达式，但这个表达式太复杂，很难手动计算导数的解析形式。或者目标函数是随时变化的，无法事先求出。一个典型的例子就是神经网络，每次我们更改网络的层数、维度、激活函数、损失函数等等，都会使得目标函数发生变化。手动求导是不可行的，数值求导虽然可行但精度难以保证，且巨大的自变量维度会导致雅可比矩阵计算量突破天际。于是，自动求导的出现解决了这个问题。

很多人搞不清楚数值求导、符号求导和自动求导的区别。本文没有提到符号求导，这是一种直接用计算机计算导数解析解的方法，存在于MATLAB、Mathematica等数值计算软件中。而自动求导不同于符号求导，它并不计算解析解，而是利用链式法则将复杂的函数拆分成一个个独立的计算单元，分别对这些小单元求导，最后再合并得到完整的导数。将函数拆分为多个计算单元后，我们称之为计算图（computational graph），这是一个有向无环图，标记了数据的流向。我们用一个例子来详细介绍这一方法。

目标函数

$f(x)=\left(x_{1} x_{2} \sin x_{3}+e^{x_{1} x_{2}}\right) / x_{3} \tag{7.26} \\$

拆分为多个计算单元

$\begin{array}{l} x_{4}=x_{1} * x_{2} \\ x_{5}=\sin x_{3} \\ x_{6}=e^{x_{4}} \\ x_{7}=x_{4} * x_{5} \\ x_{8}=x_{6}+x_{7} \\ x_{9}=x_{8} / x_{3} \end{array} \tag{7.27} \\$

按照计算顺序将其表示成有向图如下

我们最终的目的是求导数 $\nabla f(x) = \left[\begin{array}{} \frac{\partial f}{\partial x_1} & \frac{\partial f}{\partial x_2} & \frac{\partial f}{\partial x_3} \end{array}\right]$ 。以第一项 $\frac{\partial f}{\partial x_1}$ 为例，根据链式法则

$\begin{aligned} \frac{\partial f}{\partial x_1} &= \frac{\partial f}{\partial x_8} \frac{\partial x_8}{\partial x_1} \\ &= \frac{\partial f}{\partial x_8} \left( \frac{\partial x_8}{\partial x_7} \frac{\partial x_7}{\partial x_1} + \frac{\partial x_8}{\partial x_6} \frac{\partial x_6}{\partial x_1} \right) \\ &= \frac{\partial f}{\partial x_8} \left( \frac{\partial x_8}{\partial x_7} \left( \frac{\partial x_7}{\partial x_4} \frac{\partial x_4}{\partial x_1} \right) + \frac{\partial x_8}{\partial x_6} \left( \frac{\partial x_6}{\partial x_4} \frac{\partial x_4}{\partial x_1} \right) \right) \\ \end{aligned} \\$

这样分解之后，每一个需要计算的项都是相邻节点之间简单函数的导数。而且计算图从左到右对应着上式从内到外，我们只需要一次前向传播（forward sweep）即可得出结果。

自动求导就是这么简单，但还不止于此。如果我们换种方式使用链式法则，像下面这样

$\begin{aligned} \frac{\partial f}{\partial x_1} &= \frac{\partial f}{\partial x_4} \frac{\partial x_4}{\partial x_1} \\ &= \left( \frac{\partial f}{\partial x_6} \frac{\partial x_6}{\partial x_4} + \frac{\partial f}{\partial x_7} \frac{\partial x_7}{\partial x_4} \right) \frac{\partial x_4}{\partial x_1}\\ &= \left( \left( \frac{\partial f}{\partial x_8} \frac{\partial x_8}{\partial x_6} \right)\frac{\partial x_6}{\partial x_4} + \left( \frac{\partial f}{\partial x_8} \frac{\partial x_8}{\partial x_7} \right)\frac{\partial x_7}{\partial x_4} \right) \frac{\partial x_4}{\partial x_1} \\ \end{aligned} \\$

有没有发现和上面的差别？第一种链式法则是把计算图从右到左依次展开的，而第二种链式法则是把计算图从左到右依次展开的。这样的结果是，第二种方法需要先计算最右侧的导数，然后依次向左计算。因此它不仅需要一次前向传播，还需要一次反向传播（reverse sweep）。

在自动求导中，第一种方法称为前向模式（forward mode），第二种方法称为逆向模式（reverse mode）。看起来似乎逆向模式计算更复杂，而且需要保存前向传播产生的中间变量，但实际的深度学习框架中都采用的是逆向模式，这是为什么呢？

这是由目标函数的形式决定的。在深度学习以及大部分优化问题中，目标函数都是标量函数 $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$ ，自变量很多，但函数值只有一个。这种情况下，使用前向模式，需要分别对每个自变量做一次前向传播，也就是 $N$ 次前向传播。如果使用逆向模式，虽然从上面的公式看起来好像也需要 $N$ 次计算，但实际上计算的时候是从括号内部向外逐项计算的，对应到计算图上是从右向左计算。这和从左向右很不一样，因为从右向左的起始点只有一个，因此算出来的结果最终对左侧的所有输入节点都是通用的。也就是说，从右向左计算的中间节点的值可以共享给所有输入变量，从而减少计算次数。

这个道理可以推广到向量目标函数 $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^m$ 。使用前向模式，我们需要完整地计算 $N$ 次前向传播。使用逆向模式，我们需要完整地计算 $M$ 次反向传播。因此，当 $N \gg M$ 时，使用逆向模式比较好，当 $N \ll M$ 时，使用前向模式比较好。

与数值求导中的稀疏雅可比类似，自动求导时也可以考虑稀疏性。仍然是找到计算图中相互独立的输入节点，在前向传播时同时考虑这些输入即可。

总结

数值求导和自动求导是数值最优化中的基本操作，了解它们的原理，可以帮助我们更深刻地理解优化算法的性能。本文的讲解比较粗浅，旨在帮助大家了解基本概念，感兴趣的同学可以继续查阅相关资料。

见：https://zhuanlan.zhihu.com/p/109755675

免责声明！

本站转载的文章为个人学习借鉴使用，本站对版权不负任何法律责任。如果侵犯了您的隐私权益，请联系本站邮箱yoyou2525@163.com删除。

猜您在找 NDArray自动求导 python对函数值求导和求值用PyTorch自动求导 Autograd: 自动求导 Pytorch学习（一）—— 自动求导机制张量自动求导和计算图『PyTorch』第三弹_自动求导深度学习中常用的自动求导的介绍 (转)机器学习之——自动求导 TensorFlow2.0教程28：自动求导