数值求导和自动求导

2021-03-04

数值求导和自动求导

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

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

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

数值求导

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

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

使用泰勒定理展开  可得

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

进而求得导数

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

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

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

其中  。可以发现，  其实表示的是相对精度。现在回到式(7.4)，如果我们用计算机计算的结果代替理想值，即用  替代  ，用  替代  ，结果势必会变化，由式(A.58)可知

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

这意味着

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

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

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

如果我们仍按照前面的方式分别对  和  泰勒展开，可以发现式(7.7)的近似误差在  数量级。

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

稀疏雅克比

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

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

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

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

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

把有限微分公式推广到向量情况，对于雅可比矩阵中的任意一个值  ，需要取  ，其中  是第  维的单位基向量。仔细观察式(7.13)的第一列和第四列，可以发现，  只对导数中的  和  分量产生影响，  只对导数中的  、  和  分量产生影响，这说明  和  关于导数值是相互独立的。这样我们就可以取  ，一次性计算出  和  。

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

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

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

自动求导

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

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

目标函数

拆分为多个计算单元

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

我们最终的目的是求导数  。以第一项  为例，根据链式法则

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

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

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

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

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

这个道理可以推广到向量目标函数  。使用前向模式，我们需要完整地计算 次前向传播。使用逆向模式，我们需要完整地计算  次反向传播。因此，当  时，使用逆向模式比较好，当  时，使用前向模式比较好。

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

总结

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

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