Numerical Optimization Ch8. Calculating Derivatives

第八章: 优化中的导数计算

大多非线性优化和求解非线性方程的算法都需要导数的信息, 从而在计算过程中更加充分地考虑问题的性态、掌握更多关于问题的信息, 以获得较快的收敛性. 手算导数有时是可行的, 当然更加方便的做法是写一串代码让计算机帮我们算. 但在其他情形下, 函数会过于复杂, 手算显然不可取. 这时我们只能寻求自动计算或 近似计算导数的方法. 以下是较为常用的一些方法:

1. 有限差分法. 这一方法根植于Taylor定理. 通过观察在给定点$x$处基于微小扰动函数值的变动, 我们可以估计函数对于无穷小"扰动"的变化速率, 也就是导数. 例如光滑函数$f:{\mathbb{R}}^n\to \mathbb{R}$对于自变量$x$的第$i$个分量$x_i$的偏导数就可以用如下中心差分公式近似计算:$$\frac{\partial f}{\partial x_i}\approx \frac{f(x+ϵe_i)-f(x-ϵe_i)}{2ϵ},$$其中$ϵ$为一较小的正标量, $e_i$为第$i$单位向量. 利用Taylor定理, 易得这一公式计算的误差在$O({ϵ}^2)$.
2. 自动微分法. 这一方法则来源于计算机计算函数值的机理. 一般地, 计算机代码计算函数值的过程可以分解为一系列基本运算, 其中我们可以应用微积分中的链式法则. 一些软件具有计算函数值和导数值的功能. 而另一些软件则是记录计算函数在一点$x$处函数值所进行的一系列基本运算, 最后再用这些信息计算同一个点$x$处的导数值.
3. 符号微分法. 这一方法中, 计算机将使用符号运算的工具对函数$f$的代数表示进行操作, 从而产生对于每个分量导数的新的代数表示.

本章中, 我们主要讨论前两种方法——有限差分法与自动微分法.
导数的作用不仅仅在优化上. 设计优化与经济领域也经常涉及后验最佳灵敏度分析(post-optimal sensitivity analysis), 其中我们需要最优值(或最优点)对于参数或约束中微小扰动的敏感度. 导数也在非线性微分方程和模拟中发挥着重要作用.

1. 有限差分

有限差分法是一种近似计算导数的方法, 其根植于Taylor定理. 许多软件包不论用户是否能够或乐意提供计算精确导数的方法或代码, 都会自动实施有限差分. 尽管这一方法只能提供导数的近似值, 但在许多情形下得到的结果都较为令人满意.
由定义, 导数实际就是对函数对于变量上微小扰动敏感度的一种度量. 本节我们采用的方法就是在$x$的取值上加入一些微小的、有限的扰动, 再观察相应函数值的变化. 最终导数在一点处的近似值将由后者与前者的比值得到.

1.1 近似梯度

一种近似梯度向量$\nabla f(x)$的方法是, 计算$f$在$(n+1)$个点处的值, 之后进行一些基本的运算. 下面我们将介绍这种方法, 同时介绍其基于更多函数值的一种改进方案.

一种广泛使用的近似偏导$\partial f/\partial x_i$在$x$处取值的方法是前向差分法(或单边差分法), 定义为$$\frac{\partial f}{\partial x_i}(x)\approx \frac{f(x+ϵe_i)-f(x)}{ϵ}.$$而梯度就可以通过$i$遍历$1,\dots ,n$来获取. 这一过程总共需要$(n+1)$个函数值: $f(x),f(x+ϵe_1),\dots ,f(x+ϵe_n)$.
上述方法来源于Taylor定理: 当$f$二次连续可微, 我们有$$f(x+p)=f(x)+\nabla f(x{)}^{T}p+\frac{1}{2}{p}^T{\nabla }^{2}f(x+tp)p,t\in (0,1).$$若取$L$为在$x$附近$\Vert {\nabla }^{2}f(\cdot )\Vert$的上界, 则有$$|f(x+p)-f(x)-\nabla f(x{)}^{T}p|\le \frac{1}{2}L\Vert p{\Vert }^2.$$选取$p$为$ϵe_i$, 从而$\nabla f(x{)}^{T}p=ϵ\nabla f(x{)}^T{e}_i=ϵ\partial f/\partial x_i.$ 代入可得$$\left|\frac{\partial f}{\partial x_i}-\frac{f(x+ϵe_i)-f(x)}{ϵ}\right|\le {\delta }_{ϵ}=\frac{1}{2}Lϵ.$$若忽略${\delta }_{ϵ}$就有前向差分公式, 而这一误差会随着$ϵ\to 0$而趋于0.

前向差分公式应用时需要重点考量的一个方面就是参数$ϵ$的选取. 仅从误差${\delta }_{ϵ}$上看我们应当将$ϵ$取得越小越好, 但这样就忽略了计算机进行浮点运算时产生的舍入误差. 机器精度$\mathrm{u}$是计算机对两个浮点数进行单次运算时产生的相对误差的上界. 在双精度IEEE浮点运算下$\mathrm{u}\approx 1.1\times 10^{-16}$. 这些误差最终将如何影响计算的函数值将取决于$f$的计算方式, 例如一个运算公式、一个微分方程的求解器.
下面仅做粗略估计. 假设计算的$f$的相对误差以$\mathrm{u}$为上界, 从而$f(x)$和$f(x+ϵe_i)$的计算值与精确值具有如下关联:$$\begin{array}{rl}|\hat{f}(x)-f(x)|& \le \mathrm{u}{L}_f,\\ |\hat{f}(x+ϵe_i)-f(x+ϵe_i)|& \le \mathrm{u}{L}_f,\end{array}$$其中$\hat{f}(\cdot )$表示计算值, $L_f$为$|f(\cdot )|$在$x$附近的上界. 如果以计算值代入偏导数的计算, 则之前式子的误差将以$$(L/2)ϵ+2\mathrm{u}{L}_f/ϵ$$为上界. 自然地, 我们希望选取$ϵ$使得这一误差尽可能小. 易得这一对勾型函数的最小点为$${ϵ}^2=\frac{4L_f\mathrm{u}}{L}(而不是越小越好).$$若我们还假设问题是尺度良好的, 则$L_f/L$就可以用一个不太大的数控制. 因此我们可以取$$ϵ=\sqrt{\mathrm{u}}(事实上许多优化软件包均取ϵ为此值计算有限差分).$$此时, 前向差分近似的误差有上界(约为)$\sqrt{\mathrm{u}}$.

还有另外一种更加精确的近似方式——中心差分公式, 定义为$$\frac{\partial f}{\partial x_i}(x)\approx \frac{f(x+ϵe_i)-f(x-ϵe_i)}{2ϵ}.$$除去精度的优势不谈, 这种计算方式要比前向差分昂贵: 我们需要计算$f$在$(2n+1)$个点处的值. 下面我们利用Taylor定理说明其精度上的优点. 当$f$的二阶导存在且是Lipshcitz连续时, 我们有$$\begin{array}{rl}f(x+p)& =f(x)+\nabla f(x{)}^{T}p+\frac{1}{2}{p}^T{\nabla }^{2}f(x+tp)pt\in (0,1)\\ & =f(x)+\nabla f(x{)}^{T}p+\frac{1}{2}{p}^T{\nabla }^{2}f(x)p+O(\Vert p{\Vert }^3).\end{array}$$分别令$p=ϵe_i,-ϵe_i$, 则有$$f(x+ϵe_i)=f(x)+ϵ\frac{\partial f}{\partial x_i}+\frac{1}{2}{ϵ}^2\frac{{\partial }^{2}f}{\partial x_i^2}+O({ϵ}^3),$$$$f(x-ϵe_i)=f(x)-ϵ\frac{\partial f}{\partial x_i}+\frac{1}{2}{ϵ}^2\frac{{\partial }^{2}f}{\partial x_i^2}+O({ϵ}^3).$$两式相加消去$ϵ$的一次项, 整理后可得$$\frac{\partial f}{\partial x_i}(x)=\frac{f(x+ϵe_i)-f(x-ϵe_i)}{2ϵ}+O({ϵ}^2).$$我们看到这一表达式的误差在$O({ϵ}^2)$, 与前向差分的$O(ϵ)$具有很大的提升. 不过将计算$f$时的误差纳入考虑后, 类似地有最优的$ϵ$为${\mathrm{u}}^{1/3}$, 此时误差大约为${\mathrm{u}}^{2/3}$, 精度改进就不是那么明显了. 不过在某些情形下, 我们能够在计算值获取上更多精确的位数. 此时多余的耗费就变得相当有意义.

1.2 近似Jacobi矩阵

考虑向量值函数$r:{\mathbb{R}}^n\to {\mathbb{R}}^m$, 例如我们将在下一章考察的残差向量或者将在第十一章讨论的非线性方程组. 这一向量值函数的一阶导数Jacobi矩阵$J(x)$定义为$$J(x)={\left[\frac{\partial r_j}{\partial x_i}\right]}_{j=1,2,\dots ,m,i=1,2,\dots ,n}={\left[\begin{array}{c}\nabla r_1(x{)}^T\\ \nabla r_2(x{)}^T\\ ⋮\\ \nabla r_m(x{)}^T\end{array}\right]}_{m\times n},$$其中$r_j,j=1,2,\dots ,m$为$r$的分量.

1. Jacobi的某一列: 前一小节中介绍的方法自然可以直接(一次)用来计算$J(x)$的一列. 当$r$二次连续可微时, 我们可用Taylor定理得到$$\Vert r(x+p)-r(x)-J(x)p\Vert \le (L/2)\Vert p{\Vert }^2,$$其中$L$为$J$在$x$附近的Lipschitz常数.

2. Jacobi-向量乘积: 若我们需要获取Jacobi矩阵与某一给定向量$p$的乘积$J(x)p$(这在十一章使用非精确牛顿法求解非线性方程时会提到), 上述不等式立刻就给出了一种近似计算的方法:$$J(x)p\approx \frac{r(x+ϵp)-r(x)}{ϵ},$$其中误差在$O(ϵ)$. 类似地也有中心型的近似公式.

3. 完整Jacobi: 若得寸进尺, 需要求得完整的Jacobi矩阵$J(x)$, 我们可以通过一次计算一列的方式得到. 其中第$i$列的计算方法为: 令$p=ϵe_i$, $$\frac{\partial r}{\partial x_i}(x)\approx \frac{r(x+ϵe_i)-r(x)}{ϵ}.$$从上面公式可知, 要得到完整的Jacobi矩阵估计, 需要求$n+1$次向量函数值$r$. 不过当Jacobi矩阵稀疏时, 我们通常可以小得多的代价(有时甚至只需求3-4次向量函数值)得到类似的结果. 这种方法操作的关键在于, 同时估计Jacobi矩阵的多个不同列. 当然这需要对 $p$的选择加以斟酌.
   我们以一个简单的例子介绍这一方法. 考虑向量值函数$r:{\mathbb{R}}^n\to {\mathbb{R}}^n$, 定义为$$r(x)=\left[\begin{array}{c}2(x_2^3-x_1^2)\\ 3(x_2^3-x_1^2)+2(x_3^3-x_2^2)\\ 3(x_3^3-x_2^2)+2(x_4^3-x_3^2)\\ ⋮\\ 3(x_n^3-x_{n-1}^2)\end{array}\right].$$注意到$r$的每个分量均只依赖于2-3个$x$的分量, 从而可知Jacobi矩阵的每行也仅包含2-3个非零元素. 对于$n=6$的情形, Jacobi矩阵就具有如下的三对角结构:$$\left[\begin{array}{ccccccc}\times & \times & & & & & \\ \times & \times & \times & & & & \\ & \times & \times & \times & & & \\ & & \times & \times & \times & & \\ & & & \times & \times & \times & \\ & & & & \times & \times & \end{array}\right].$$现在我们要计算Jacobi的近似. 我们做一个比较.

   • 首先使用原来的"笨办法". 对$x$的第一个分量作扰动$p=ϵe_1$将只能影响$r$的第1和第2分量, 而其余分量将保持不变. 从而计算出的$\frac{\partial r}{\partial x_1}$在第3,4,5,6分量上都是0. 如果我们提前已经知道了这些分量是什么, 再去重复计算显然是不应该的.
   • 再使用"新方法". 相比于$ϵe_1$这样的只能改变2个分量的扰动向量, 我们更想要寻求能够影响后4个分量的扰动向量. 这将成为我们同时扰动多列方法的基础. 在这个例子中不难看出$ϵe_4$就(差不多)具有这样的性质: 它会改变$r$中第3,4,5分量的值, 但会保持第1,2分量的值不变. 进一步地, 我们可以说$ϵe_1$和$ϵe_4$带来的扰动是互补影响的.
     为严格地进行数学推导, 我们令$$p=ϵ(e_1+e_4)$$并注意到$$r(x+p{)}_{1,2}=r(x+ϵ(e_1+e_4){)}_{1,2}=r(x+ϵe_1{)}_{1,2},$$$$r(x+p{)}_{3,4,5}=r(x+ϵ(e_1+e_4){)}_{3,4,5}=r(x+ϵe_4{)}_{3,4,5}.$$代入可得$$r(x+p{)}_{1,2}=r(x{)}_{1,2}+ϵ[J(x)e_1{]}_{1,2}+O({ϵ}^2).$$重新整理式子, 我们就可计算Jacobi矩阵的(1,1)元和(2,1)元:$$\left[\begin{array}{c}\frac{\partial r_1}{\partial x_1}(x)\\ \frac{\partial r_2}{\partial x_1}(x)\end{array}\right]=[J(x)e_1{]}_{1,2}\approx \frac{r(x+p{)}_{1,2}-r(x{)}_{1,2}}{ϵ}.$$类似地,$$\left[\begin{array}{c}\frac{\partial r_3}{\partial x_4}(x)\\ \frac{\partial r_4}{\partial x_4}(x)\\ \frac{\partial r_5}{\partial x_4}(x)\end{array}\right]=[J(x)e_4{]}_{3,4,5}\approx \frac{r(x+p{)}_{3,4,5}-r(x{)}_{3,4,5}}{ϵ}.$$这就是说, 我们可通过只扰动一次、只计算一次函数值$r(x+ϵ(e_1+e_4))$便可同时得到Jacobi矩阵两列的估计.

   我们也可以经济地估计$J(x)$剩下的几列. 比如第2,5列就可通过设置$p=ϵ(e_2+e_5)$同时得到, 第3,6列可通过设置$p=ϵ(e_3+e_6)$同时得到. 总之, 我们仅需计算3次函数值.
   事实上, 对于上述例子的任意$n$的情形, 计算3次函数值都足够估计整个Jacobi了. 对应的扰动向量$p$为$$\begin{gathered} p=ϵ(e_1+e_4+e_7+e_{10}+\cdots ), \\ p=ϵ(e_2+e_5+e_8+e_{11}+\cdots ), \\ p=ϵ(e_3+e_6+e_9+e_{12}+\cdots ). \end{gathered}$$这种设置扰动向量的方式也可启发其他向量值函数的情形.
   事实上, 选取扰动向量的算法可以用图和图着色的语言简洁的表示出来. 瑞昱任一向量值函数$r:{\mathbb{R}}^n\to {\mathbb{R}}^m$, 我们可建立有$n$个顶点的列关联图$G$, 其中若$r$有分量同时依赖于$x_i,x_k$(即$J(x)$的第$i,k$列均在某个第$j$行有非零元), 对应地在$G$中节点$i,k$之间就有边连接. 前面的$n=6$的例子对应的列关联图如下所示.
   

   下面以某种规则对$G$进行着色, 规则: 任意两相邻接的顶点颜色不同. 最终我们将根据颜色选取扰动向量: 若节点$i_1,i_2,\dots ,i_l$有相同的颜色, 则就有扰动向量$p=ϵ(e_{i_1}+e_{i_2}+\cdots +e_{i_l})$.
   一般地, 满足上述规则的图着色方案不止一种. 最简单的方式就是对每个节点均赋以不同的颜色, 不过这样就对应了$n$个扰动向量, 这样的效率也是最低的. 找寻一个图所需的最少着色数是一个NP-完全问题¹, 不过我们可以使用低成本的算法求得近似最优的解. Newsam和Ramsdell²表示通过考虑更广泛的一类扰动向量, 我们甚至可能仅用不超过$n_z$个函数值便可得到完整Jacobi矩阵的近似. 这里$n_z$是$J(x)$中最大的行非零元个数.
   对于一些具有良好结构(已经被深入研究了)的Jacobi矩阵(大量可见于微分算子的离散或能够得到带型Jacobi的情形), 最优的着色方案是已知的. 例如对于上图中的三对角结构, 3种颜色的方案就是最佳方案.

1.3 近似Hessian矩阵

有时, 我们可能可以计算梯度$\nabla f(x)$而不能计算Hessian矩阵${\nabla }^{2}f(x)$. 我们当然可以通过上面介绍的方法计算(将$r$换成$\nabla f$). 通过图着色, 我们可以以极小的代价近似稀疏的Hessian. 然而这些做法都忽略了Hessian矩阵对称性的特征, 且这一忽略可能最终导致算得的Hessian并不对称. 我们当然可以通过加转置除2的方式获取Hermite部分. 但这种方式所带来的影响难以估计. 下面我们介绍一种将${\nabla }^{2}f(x)$的对称性纳入考量的差分法.

1. Hessian-向量乘积: 一些重要的算法——比如第七章描述的Newton-CG算法——并不需要完整的Hessian矩阵, 而仅需要我们提供矩阵-向量乘积${\nabla }^{2}f(x)p,\forall p$. 我们可以用Taylor定理再一次地估计这一乘积: 当$f$的二阶导数存在且在$x$附近Lipschitz连续时, 我们有$$\nabla f(x+ϵp)=\nabla f(x)+ϵ{\nabla }^{2}f(x)p+O({ϵ}^2),$$从而$${\nabla }^{2}f(x)p\approx \frac{\nabla f(x+ϵp)-\nabla f(x)}{ϵ}.$$这一表达式的误差在$O(ϵ)$, 此时额外的代价是计算在$x+ϵp$处的梯度. 上面这一公式对应之前的前向差分法. 自然若计算$\nabla f(x-ϵp)$我们也有中心差分法.
2. 完整Hessian: 当梯度都无法计算时, 我们只能只用函数值来近似Hessian矩阵的元素. 此时依然要用到Taylor定理. 此时令$p=ϵe_i,ϵe_j,ϵ(e_i+e_j)$, 则$$\begin{gathered} f(x+ϵe_i)=f(x)+ϵ\nabla f(x{)}^T{e}_i+{ϵ}^2\frac{1}{2}{e}_i^T{\nabla }^{2}f(x)e_i+O({ϵ}^3), \\ f(x+ϵe_j)=f(x)+ϵ\nabla f(x{)}^T{e}_j+{ϵ}^2\frac{1}{2}{e}_j^T{\nabla }^{2}f(x)e_j+O({ϵ}^3), \\ f(x+ϵ(e_i+e_j))=f(x)+ϵ\nabla f(x{)}^T(e_i+e_j)+{ϵ}^2\frac{1}{2}(e_i+e_j{)}^T{\nabla }^{2}f(x)(e_i+e_j)+O({ϵ}^3). \end{gathered}$$从而经整理可得$$\frac{{\partial }^{2}f}{\partial x_i\partial x_j}(x)=\frac{f(x+ϵ(e_i+e_j))-f(x+ϵe_i)-f(x+ϵe_j)+f(x)}{{ϵ}^2}+O(ϵ).$$若我们想用这个公式近似Hessian的每一个元素, 我们就需要计算$f$在$x+ϵ(e_i+e_j),\forall i,j$处的取值(至多$n(n+1)/2$个不同取值), 以及另外$n$个点$x+ϵe_i,i=1,2,\dots ,n$处的值. 如果Hessian稀疏, 且我们知道哪些元素是零的话, 当然这一步就可以跳掉.

1.4 近似稀疏Hessian矩阵

之前我们提到, 可以对$\nabla f$进行差分以得到Hessian的近似. 下面我们介绍如何在Hessian稀疏时,利用Hessian的对称性来减少获取完整近似所需的扰动向量$p$的个数. 关键的一点就是, 对于$[{\nabla }^{2}f(x){]}_{i,j}$的估计同样也是对$[{\nabla }^{2}f(x){]}_{j,i}$的估计.
我们以一个简单函数$f:{\mathbb{R}}^n\to \mathbb{R}$说明, 其定义为$$f(x)=x_1\sum _{i=1}^n{i}^2{x}_i^2.$$利用微积分, 我们知道${\nabla }^{2}f$具有"爪型"结构, 以$n=6$为例:$$\left[\begin{array}{cccccc}\times & \times & \times & \times & \times & \times \\ \times & \times & & & & \\ \times & & \times & & & \\ \times & & & \times & & \\ \times & & & & \times & \\ \times & & & & & \times \end{array}\right].$$若建立对$\nabla f$的列关联图, 则我们会发现图中每两个节点都是相邻接的, 这是因为${\nabla }^{2}f$的第一行无0元. 由图着色的规则, 我们就需要给不同的节点赋不同的颜色, 即需要计算$\nabla f$在$n+1$个点处的值.

我们可以利用对称性建立更加高效的方法. 假设我们先令$p=ϵe_1$估计${\nabla }^{2}f(x)$的第一列. 由对称性, 同样的估计也适用于${\nabla }^{2}f(x)$的第一行. 此时我们只剩下${\nabla }^{2}f(x{)}_{22},{\nabla }^{2}f(x{)}_{33},\dots ,{\nabla }^{2}f(x{)}_{66}$没有估计. 但注意到此时剩下的节点2-6构成的子图中, 节点之间互不邻接. 这就是说, 我们可以赋给它们相同的颜色, 即令扰动向量为$$p=ϵ(e_2+e_3+\cdots +e_6)=ϵ(0,1,1,1,1,1{)}^T.$$注意到$\nabla f$的第2个分量与未知向量的第3,4,5,6个分量无关, 而第3个分量与未知向量的第2,4,5,6个分量无关, 以此类推. 因此, 对于第$i$个分量, $$\nabla f(x+p{)}_i=\nabla f(x+ϵ(e_2+e_3+\cdots +e_6){)}_i=\nabla f(x+ϵe_i{)}_i.$$对这些单独的分量应用前向差分公式可得$$\frac{{\partial }^{2}f}{\partial x_i^2}(x)\approx \frac{\nabla f(x+ϵe_i{)}_i-\nabla f(x{)}_i}{ϵ}=\frac{\nabla f(x+ϵp{)}_i-\nabla f(x{)}_i}{ϵ},i=2,3,\dots ,6.$$这样, 我们利用对称性, 仅需要计算$\nabla f$在$x$和另外两个点处的值就可以估计出整个稀疏的Hessian.

图着色技巧可以再次用来经济地选取我们的扰动向量$p$. 此处, 我们使用邻接图代替之前使用的列关联图. 邻接图具有$n$个顶点, 且当$i̸=k$, ${\partial }^{2}f(x)/(\partial x_i\partial x_k)̸=0$时在节点$i,k$邻接. 此时的图着色规则也比之前的更加复杂: 我们不仅要求邻接的节点具有不同的颜色, 且要求任一图中长度为3的路径包含至少3种颜色. 换句话说, 若在图中有节点$i_1,i_2,i_3,i_4$且有边$(i_1,i_2),(i_2,i_3),(i_3,i_4)$, 则在给这4个节点着色时必须至少用3种不同的颜色. 扰动向量的构造与之前相同: 若节点$i_1,i_2,\dots ,i_l$具有相同颜色, 则设置扰动向量为$$p=ϵ(e_{i_1}+e_{i_2}+\cdots +e_{i_l}).$$

2. 自动微分

自动微分是一类技术的统称, 这类技术使用函数的计算表示来产生导数的解析值. 其中一些直接通过对函数值进行操作获取导数. 其他的则会保留计算某点$x$处函数值的计算过程并回过头来利用这些信息计算$x$处的导数值.

自动微分的详细理论、软件工具和应用可见 www.autodiff.org.

自动微分技术建立在这样的观察上: 任何函数, 不论多复杂, 想要计算它的值都是通过一系列的一元或二元运算. 二元运算包括加法、乘法、除法和求幂. 一元运算有比如三角函数、指数函数和对数函数(作为媒介). 自动微分的另一个共同的要素就是, 它们均会使用微积分中的链式法则.
有两种基本的自动微分的模式: 前向模式与反向模式. 它们的区别可以用一个简单的例子说明. 我们还将利用这个例子阐述这些技术如何延展到一般函数上, 包括向量值函数.

2.1 一个例子

考虑以下具有3个变量的函数:$$f(x)=(x_1{x}_2\sin x_3+e^{x_1{x}_2})/x_3.$$下图表示了计算这样的函数值可以分解为怎样的基本运算, 并也表明了这些运算的偏序如何.

例如, 乘法$x_1\ast x_2$必须要先于指数运算$e^{x_1{x}_2}$, 否则我们会得到错误的结果$e^{x_1}{x}_2$. 上图中引入了中间变量$x_4,x_5,\dots$来储存计算过程中的中间值; 它们与图中最左边的独立变量$x_1,x_2,x_3$不同. 我们可以算术形式表示计算$f$值的过程:$$\begin{array}{rl}{x}_4& =x_1\ast x_2,\\ x_5& =\sin x_3,\\ x_6& =e^{x_4},\\ x_7& =x_4\ast x_5,\\ x_8& =x_6+x_7,\\ x_9& =x_8/x_3.\end{array}$$上图中最后的节点$x_9$则储存了函数值$f(x)$. 以图论的语言说, 当有一有向弧从节点$i$指向节点$j$, 我们称节点$i$为节点$j$的父节点, 节点$j$为节点$i$的子节点. 任一节点只有在它所有的父节点都已知的条件下才可以计算, 因此计算过程在图中从左向右. 这样的计算过程称作是前推(forward sweep). 需要指明的是, 用于做自动微分的软件工具并不需要用户将计算函数值的代码分解成计算一个又一个中间值的小组分. 中间值的体量识别和构建都可靠软件工具本身显式或隐式地完成.

2.2 前向模式

在自动微分的前向模式中, 我们将推进和计算每个中间变量$x_i$对某一给定方向$p\in {\mathbb{R}}^n$的方向导数, 同时计算$x_i$本身. 对于上述3个独立变量的例子, 我们使用下述表达式来表示每个变量对$p$的方向导数:$$D_p{x}_i\stackrel{def}{=}(\nabla x_i{)}^{T}p=\sum _{j=1}^3\frac{\partial x_i}{\partial x_j}{p}_j,i=1,2,\dots ,9,$$其中$\nabla$表示对3个独立变量求梯度. 我们的最终目标就是计算${\nabla }_p{x}_9$, 即$\nabla f(x{)}^{T}p$. 注意到独立变量$x_i,i=1,2,3$的方向导数分别为$p_i,i=1,2,3$. 方向$p$称为种子向量(seed vector).
任意节点处只要得到了节点处的值$x_i$, 我们就可用链式法则计算$D_p{x}_i$. 例如我们已知$x_4,D_p{x}_4,x_5,D_p{x}_5$, 而由图我们知道$x_7=x_4{x}_5$. 由链式法则, $$\nabla x_7=\frac{\partial x_7}{\partial x_4}\nabla x_4+\frac{\partial x_7}{\partial x_5}\nabla x_5=x_5\nabla x_4+x_4\nabla x_5.$$上式两边与$p$做内积可得$$D_p{x}_7=\frac{\partial x_7}{\partial x_4}{D}_p{x}_4+\frac{\partial x_7}{\partial x_5}{D}_p{x}_5=x_5{D}_p{x}_4+x_4{D}_p{x}_5.$$于是方向导数$D_p{x}_i$就被逐个获取, 并最终得到$$D_p{x}_9=D_{p}f=\nabla f(x{)}^{T}p.$$
前向模式的原理是很直接的, 但它的实际操作和计算需求又如何呢?

1. 首先, 我们再次强调我们在实际应用时不需要构建计算图、将计算过程组分化或者是定义中间变量. 自动微分的软件工具会隐式、自动地帮助我们实施. 我们也不需要一直储存计算图中每一个节点的信息$x_i,D_p{x}_i$. 这是因为当一个节点的所有子节点全都计算完毕后, 它的$x_i,D_p{x}_i$也就不再需要了.
2. 实际操作中的关键是逐步逐个地计算$x_i,D_p{x}_i$. 自动微分的软件会在计算过程中对任何标量$w$关联另一个标量$D_{p}w$. 只要$w$在计算过程中有使用, 软件就会自动对这一计算过程附带关联上梯度向量$D_{p}w$. 例如, 若$w$被用在一个除法运算中, $$z\leftarrow \frac{\partial w}{\partial y},$$此时就会调用$w,z,D_{p}w,D_{p}y$计算方向导数$D_{p}z$:$$D_{p}z\leftarrow \frac{1}{y}{D}_{p}w-\frac{w}{y^2}{D}_{p}y.$$为获取完整的梯度向量, 我们可以同时对$n$个种子向量$p=e_1,e_2,\dots ,e_n$进行操作. 此时$p=e_j$对应$D_{p}f=\partial f/\partial x_j,j=1,2,\dots ,n$. 而从上面的例子我们发现, 计算$f$和$\nabla f$的额外耗费可能会相当大. 就此例而言, 单单一个$w$除以$y$就需要我们在计算$D_{e_j}z$时计算$2n$次乘法和$n$次加法,$j=1,2,\dots ,n$. 而这仅仅是一个除法, 对于实际更复杂的关系将会使计算的耗费难以估计(这里还有获取调用和储存数据的耗费). 储存量可能会随着与$n$同阶的因子增长, 这是因为我们现在需要对中间变量$x_i$存储$n$个额外的标量$D_{e_j}{x}_i,j=1,2,\dots ,n$. 此时若能了解到这些量有很多(甚至哪些)是零的话, 可能就能节省一些存储空间. 这点在计算初期(也就是计算图的左端)尤为重要, 此时我们就可以用稀疏的数据结构储存向量$D_{e_j}{x}_i,j=1,2,\dots ,n$.

自动微分前向模式可以一种预编译的方式实施. 它会将函数值计算的代码扩展为也能计算导数的代码. 另一种方式是使用C++等语言中可用的运算符重载工具, 以上述方式透明地扩展数据结构和运算.

2.3 反向模式

自动微分的反向模式并不同时计算函数值和梯度, 而是在完成函数值的计算后, 再反过来看$f$对各个变量$x_i$——诸如独立变量和中间变量——的偏导数. 这一步是通过对计算图的反推(reverse sweep)实现的. 在这一过程的最后, 梯度向量$\nabla f$就从偏导数$\partial f/\partial x_i,i=1,2,\dots ,n$组装起来.
不像前向模式中使用$D_p{x}_i$, 反向模式将对图中每个节点关联一个标量变量${\bar{x}}_i$; 在反推过程中, 关于偏导数$\partial f/\partial x_i$的信息就储存在${\bar{x}}_i$中. ${\bar{x}}_i$有时称作伴随变量(adjoint variables), 其初始值设成0, 除了图最右边的节点$N$我们设${\bar{x}}_N=1$. 这一设置是合理的, 这是因为$x_N$储存了最终的函数值$f$, 从而$\partial f/\partial x_N=1$.
反推的过程也是基于链式法则: 对任一节点$i$, 偏导$\partial f/\partial x_i$可以从其子节点$j$的偏导$\partial f/\partial x_j$中计算:$$\frac{\partial f}{\partial x_i}=\sum _{j\mathrm{a}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{o}\mathrm{f}i}\frac{\partial f}{\partial x_j}\frac{\partial x_j}{\partial x_i}.$$对每个节点$i$, 只要上式右端有一项已知, 我们就把它加到${\bar{x}}_i$上; 即$${\bar{x}}_i+=\frac{\partial f}{\partial x_j}\frac{\partial x_j}{\partial x_i}.$$当节点$i$所有的子节点全都算完, 我们就有${\bar{x}}_i=\partial f/\partial x_i$, 此时我们声明节点$i$已完成计算(finalized). 之后, 节点$i$就将作为其父节点的子节点进一步计算. 这一过程以这种方式进行直至所有的节点均以完成计算.

在反推的过程中, 我们的计算仅涉及数值, 而不涉及公式或者计算变量$x_i$、偏导$\partial f/\partial x_i$的计算机代码. 而在前推(注意反向模式也有前推)——计算$f$——我们不仅需要计算每个变量$x_i$, 还需要计算并储存每个偏导$\partial x_j/\partial x_i$的数值. 每个偏导都与计算图中一段特定的弧相关. 在反推过程中就要用到在前推过程中计算的数值$\partial x_j/\partial x_i$.
我们以之前2.1中的例子再一次地阐释反向模式. 在下图中我们以一个特定点$x=(1,2,\pi /2{)}^T$填充. 每个节点上的数值代表了之间变量的计算值, 有向边上的数值则代表了偏导.

设置${\bar{x}}_i=0,{\bar{x}}_9=1$. 此时节点9已经完成计算(由于没有子节点).
节点9是节点3和节点8的子节点, 从而更新${\bar{x}}_3,{\bar{x}}_8$的值:$$\begin{array}{rl}{\bar{x}}_3& +=\frac{\partial f}{\partial x_9}\frac{\partial x_9}{\partial x_3}=-\frac{2+e^2}{(\pi /2{)}^2}=\frac{-8-4e^2}{{\pi }^2},\\ {\bar{x}}_8& +=\frac{\partial f}{\partial x_9}\partial x_9\partial x_8=\frac{1}{\pi /2}=\frac{2}{\pi }.\end{array}$$此时节点3还未完成计算, 而节点8已经完成, 从而$\frac{\partial f}{\partial x_8}=2/\pi$. 接着更新节点8的两个父节点:$$\begin{array}{rl}{\bar{x}}_6& +=\frac{\partial f}{\partial x_8}\frac{\partial x_8}{\partial x_6}=\frac{2}{\pi },\\ {\bar{x}}_7& +=\frac{\partial f}{\partial x_8}\frac{\partial x_8}{\partial x_7}=\frac{2}{\pi }.\end{array}$$此时节点6,7已完成计算. 之后更新节点4,5. 最终, 所有节点均计算完成, 在节点1,2,3上有$$\left[\begin{array}{c}{\bar{x}}_1\\ {\bar{x}}_2\\ {\bar{x}}_3\end{array}\right]=\nabla f(x)=\left[\begin{array}{c}(4+4e^2)/\pi \\ (2+2e^2)/\pi \\ (-8-4e^2)/{\pi }^2\end{array}\right],$$从而导数计算完成.

1. 反向模式的主要优点在于, 其对于标量函数$f:{\mathbb{R}}^n\to \mathbb{R}$的计算复杂度较低, 而计算梯度的额外计算量至多是单独计算函数值的计算量的4-5倍. 以上例中更新节点3,8的过程为例, ${\bar{x}}_3$的更新用了2次乘法、1次除法和一次加法, 而${\bar{x}}_8$的更新用了1次除法和1次加法. 这大概就是在前推过程计算函数值时的1次除法运算量的五倍.
   2.2中提到, 前向模式可能需要多于$n$倍计算函数值的计算量用于计算梯度$\nabla f$, 这使得它相比于反向模式毫无竞争力. 不过我们在下节会提到, 当我们考虑向量值函数$r:{\mathbb{R}}^n\to {\mathbb{R}}^m$时, 前向模式和反向模式的相对计算耗费随着$m$的增大将逐渐相近.
2. 反向模式也有一个很显然的缺陷, 那就是它需要储存整个计算图用于反推过程. 原则上, 计算图的存储并不难以实施: 每当在前推进行一次基本运算, 我们就可以构造储存一个新的节点用于包容中间的计算结果, 并将这一节点指向(1个或2个)父节点, 同时计算与这些边相关的偏导数. 在反推过程中, 节点就可以写入的逆序读取. 构造和写入计算图的过程可看做基本运算(通过运算符重载)的直接扩展. 反推过程或者是梯度计算的过程可直接以函数形式调用.
   不幸的是, 储存计算图所需计算量可能甚是庞大. 如果每个节点都用20个字节储存, 在每秒1亿次浮点运算的计算机上需要1秒评估时间的函数可最多产生2千兆字节大小的图. 若加以额外的计算(如在计算图的一部分上进行部分前推和后推), 则存储量就可减少.

2.4 向量值函数与部分可分性

我们已经讨论了对于一般标量值函数$f:{\mathbb{R}}^n\to \mathbb{R}$的自动微分方法. 而在(第十章就将接触的)非线性最小二乘问题和(第十一章中的)非线性方程组中, 我们必须处理向量值函数$r:{\mathbb{R}}^n\to {\mathbb{R}}^m$, 其中有$m$个分量$r_j,j=1,2,\dots ,m$. 此时计算图的最右列就有$m$个节点, 它们都没有任何的子节点, 这与之前标量值函数对应计算图中只有1个节点形成鲜明对比. 显然前向模式和反向模式可以直接推广至近似计算Jacobi矩阵$J(x)$.

向量值函数的自动微分除了在最小二乘与非线性方程问题中有应用之外, 也是处理部分可分函数的有效方法. 根据上一章最后部分介绍的, 部分可分性常见于大规模优化, 而我们可以巧妙地利用这一结构应用高效的拟牛顿方法. 当前对于给定函数$f$, 我们能够使用一些自动化的工具探测它的部分可分表示. 这为我们挖掘这一性质具有的高效性提供了可能. 我们不需要再向用户索要更多关于函数的信息.
在最简单的情形下, 函数$f$部分可分, 如果我们能把它表示成以下形式$$f(x)=\sum _{i=1}^{ne}{f}_i(x),$$其中每个基本函数$f_i(\cdot )$仅依赖于$x$的少数分量. 若我们由部分可分性构建向量值函数$r$, 即$$r(x)=\left[\begin{array}{c}{f}_1(x)\\ f_2(x)\\ ⋮\\ f_{ne}(x)\end{array}\right],$$则$$\nabla f(x)=J(x{)}^{T}e,$$其中$e=(1,1,\dots ,1{)}^T$. 由部分可分性, $J(x)$的大多列只有少数非零元. 这样我们就能通过下一小节中的图着色方法高效的计算$J(x)$. 而$\nabla f(x)$那可以用上一公式获得.
在约束优化中, 同时计算目标函数$f$和约束函数$c_i,i\in \mathcal{I}\cup \mathcal{E}$往往能带来便利. 如此一来, 我们还能够有效利用相同的表达式(即在计算图中共享相同的中间节点)有效减少总的负载(亦如之前例子中, 节点4由节点6,7所共享). 此时向量值函数$r$定义为$$r(x)=\left[\begin{array}{c}f(x)\\ [c_j(x){]}_{j\in \mathcal{I}\cup \mathcal{E}}\end{array}\right].$$

2.5 计算向量值函数的Jacobi矩阵

1. 计算Jacobi-向量乘积:

   1. 应用前向模式时, 向量值函数与标量函数是相同的. 给定种子向量$p$, 我们持续地将$D_p{x}_i$关联到储存中间变量$x_i$的节点上. 最右端的节点则会储存$D_p{r}_j=(\nabla r_j{)}^{T}p,j=1,2,\dots ,m$. 将这$m$个量组装起来, 就得到了$J(x)p$. 对于$m=1$的情形(即标量值函数), 我们可设置$p=e_1,e_2,\dots ,e_n$并同时计算$n$个量$D_{e_j}{x}_i$得到Jacobi矩阵的估计. 而对于稀疏Jacobi, 我们可采用之前有限差分法中的图着色技巧来恰当选取种子向量$p$. 相较于$r$的单次函数值计算, 计算量上的增加的倍数大约是使用的种子向量的个数.
   2. 对向量值函数应用反向模式的关键在于种子向量$q\in {\mathbb{R}}^m$的选取, 并将反向模式应用于标量值函数$r(x{)}^{T}q$. 这一过程的结果是$$\nabla [r(x{)}^{T}q]=\nabla \left[\sum _{j=1}^m{q}_j{r}_j(x)\right]=J(x{)}^{T}q.$$与前向模式获取Jacobi-向量乘积不同, 反向模式会给出Jacobi转置-向量乘积. 此时需设置储存$r_1,r_2,\dots ,r_m$的$m$个依赖节点中的变量${\bar{x}}_i$为$q$的$m$个分量. 反推到最后, 储存独立变量$x_1,x_2,\dots ,x_n$的节点中就有$$\frac{d}{d{x}_i}[r(x{)}^{T}q],i=1,2,\dots ,n,$$这些就是$J(x{)}^{T}q$的分量.

2. 计算完整Jacobi: 与之前相同, 我们通过设置$m$个单位向量$q=e_1,e_2,\dots ,e_m$获取完整的一般Jacobi估计. 而对于稀疏的Jacobi, 我们可应用之前的图着色技巧寻找更少的种子向量——唯一不同在于, 此时的图和着色策略均基于$J(x{)}^T$, 而不是$J(x)$. 同时, 计算量的增加倍数也将超过5(注意5是对于标量函数的通常的上界). 储存计算图所需的空间与标量情形差不多. 我们仅需储存图的拓扑信息以及每条边上附带的偏导数.

   前向和反向模式可以混合使用以计算$J(x)$的所有元素. 我们可为前向模式选取种子向量$p$的集合得到$J$的某些列, 再为反向模式选取种子向量$q$的集合得到$J$的包含剩下的元素的行. 也就是说, 我们可以通过行列分割的方式将前向与反向有机地结合在一起.

最后我们说明, 对于某些算法, 我们不需要提供Jacobi$J(x)$的完整信息. 例如对于非线性方程组使用的非精确牛顿法只需要提供对于任一向量$p$Jacobi与向量的乘积$J(x)p$. 而这一点用一次前推就能做到, 计算量与单独计算函数值差不多.

2.6 计算Hessian矩阵: 前向模式

至今, 我们已经介绍了如何使用前向模式与反向模式计算标量值和向量值函数的一阶导数. 下面我们概述如何扩展这些方法计算标量函数$f$的Hessian矩阵${\nabla }^{2}f$, 以及计算Hessian与给定向量$p$的Hessian-向量乘积${\nabla }^{2}f(x)p$.
回忆在前向模式中, 我们用到了$D_p{x}_i$. 现对于给定种子向量对$p,q\in {\mathbb{R}}^n$, 我们对计算图中的节点$i$定义另一个标量$$D_{pq}{x}_i=p^T({\nabla }^2{x}_i)q.$$我们可用前推基于$x_i,D_p{x}_i$计算这些量. 在储存独立变量$x_i$的节点上初始值$D_{pq}$为0(这是合理的). 当前推完成, 计算图中最右端的值$D_{pq}{x}_i$就是$p^T{\nabla }^{2}f(x)q$.
计算$D_{pq}{x}_i$的公式同样来源于链式法则. 例如, 若$x_i$是通过其两个父节点相加得来, $x_i=x_j+x_k$, 则在$D_p{x}_i$和$D_{pq}{x}_i$上对应的加和运算为:$$D_p{x}_i=D_p{x}_j+D_p{x}_k,D_{pq}{x}_i=D_{pq}{x}_j+D_{pq}{x}_k.$$其他的二元运算$-,\times ,÷$也可类似处理. 若$x_j$是通过对$x_j$做酉变换得到的, 则我们有$$\begin{array}{rl}{x}_i& =L(x_j),\\ D_p{x}_i& =L^{\prime }(x_j)(D_p{x}_j),\\ D_{pq}{x}_i& =L^{\prime \prime }(x_j)(D_p{x}_j)(D_q{x}_j)+L^{\prime }(x_j)D_{pq}{x}_j.\end{array}$$我们可以看到第三个式子中$D_{pq}{x}_i$依赖于一阶导$D_p{x}_i,D_q{x}_i$, 因此这些量也必须在前推的时候累积起来.

1. 完整Hessian: 若想计算一般完整的Hessian, 我们需要种子向量对$(p,q)$取遍所有可能的单位向量$(e_j,e_k),j=1,2,\dots ,n,k=1,2,\dots ,j$的组合, 共$n(n+1)/2$对. 当我们知晓${\nabla }^{2}f(x)$的稀疏结构后, 我们就只需要对${\nabla }^{2}f(x)$可能非零的那些位置$(j,k)$计算$D_{e_j{e}_k}{x}_i$了.
   总的计算量大约是单独计算函数值计算量的$c(1+n+N_z({\nabla }^{2}f))$倍, 其中$N_z({\nabla }^{2}f)$为${\nabla }^{2}f$中我们需要计算的元素个数. 这一数目反映了对于$N_z({\nabla }^{2}f)$个向量对$(e_j,e_k)$需要计算的量($x_i,D_{e_j}{x}_i(j=1,2,\dots ,n),D_{e_j{e}_k}{x}_i$)有多少. 小量$c$的存在则是因为更新$D_p{x}_i,D_{pq}{x}_i$需要比更新$x_i$耗费更多的计算量(比如上面酉变换的例子). 至于存储量, 图中每个节点均需要存储$1+n+N_z({\nabla }^{2}f)$个量, 不过注意当一个节点的所有节点都计算完后它的存储是可以释放的.
2. Hessian-向量乘积: 当我们不需要完整的Hessian, 而仅需矩阵-向量乘积(例如第七章的Newton-CG)时, 总的计算量必然会减小. 给定向量$q\in {\mathbb{R}}^n$, 我们使用上面的方法计算一阶导数$D_{e_1}{x}_i,\dots ,D_{e_n}{x}_i$和$D_q{x}_i$以及二阶导数$D_{e_{1}q}{x}_i,\dots ,D_{e_{n}q}{x}_i$. 于是最后的节点会依次储存$$e_j^T({\nabla }^{2}f(x))q=[{\nabla }^{2}f(x)q{]}_j,j=1,2,\dots ,n,$$而这就是${\nabla }^{2}f(x)q$的分量. 由于在前推时除了$x_i$还需要另外计算$2n+1$个量, 因此在计算量上将仅增加$2n$倍.
3. 稀疏Hessian: 计算稀疏Hessian的一种替代方案则基于单变量函数一阶和二阶导数的前向模式传播. 为说明, 注意Hessian的第$(i,j)$元可以表示成$$\begin{array}{rl}[{\nabla }^{2}f(x){]}_{ij}& =e_i^T{\nabla }^{2}f(x)e_j\\ & =\frac{1}{2}[(e_i+e_j{)}^T{\nabla }^{2}f(x)(e_i+e_j)-e_i^T{\nabla }^{2}f(x)e_i-e_j^T{\nabla }^{2}f(x)e_j].\end{array}$$当所有的节点$x_k$处的二阶导数$D_{pp}{x}_k,p=e_i,e_j,e_i+e_j$均在前推时计算完成时, 我们就可以用上面的插值公式计算$[{\nabla }^{2}f(x){]}_{ij}$.
   这种方法的一大好处就是, 我们不需要再考虑计算$D_{pq}{x}_k,p̸=q$了. 而每个$D_{pp}{x}_k$都是$x_l,D_p{x}_l,D_{pp}{x}_l$的函数, 其中节点$l$为$k$的父节点.
   注意, 若我们定义单变量函数$\psi$为$$\psi (t)=f(x+tp),$$则$D_{p}f,D_{pp}f$的值就是$\psi$的一阶和二阶导数在$t=0$处的值, 即$$D_{p}f=p^T\nabla f(x)={\psi }^{\prime }(t){|}_{t=0},D_{pp}f=p^T{\nabla }^{2}f(x)p={\psi }^{\prime \prime }(t){|}_{t=0}.$$将这一方法推广至计算三阶、四阶乃至更高阶的导数也是有可能的. 类似的插值公式可以与$psi$的更高阶的导数混合使用, 其中需恰当选取向量$p$(由一些单位向量的和组成).

2.7 计算Hessian矩阵: 反向模式

我们也可以设计基于反向模式的算法计算Hessian-向量乘积${\nabla }^{2}f(x)q$或者是整个Hessian矩阵${\nabla }^{2}f(x)$. 一种计算${\nabla }^{2}f(x)q$的算法如下: 我们首先使用前向模式通过前推时累加两个变量$x_i,D_q{x}_i$, 计算$f$和$\nabla f(x{)}^{T}q$. 之后以反向模式应用于$\nabla f(x{)}^{T}q$. 反推到最后, 计算图中代表独立变量的的节点$i=1,2,\dots ,n$就会含有$$\frac{\partial }{\partial x_i}(\nabla f(x{)}^{T}q)=[{\nabla }^{2}f(x)q{]}_i,i=1,2,\dots ,n.$$以这种方式计算${\nabla }^{2}f(x)q$的计算量相较于计算$f$增加的并不多, 且独立于$n$. 通常, 前向模式中计算$f$和$\nabla f(x{)}^{T}q$所需的计算量为单独计算$f$的一个小倍数, 而反向模式则至多引入倍数因子5. 而总的增加倍数大约是12. 若我们需要整个${\nabla }^{2}f(x)$, 就得对$q=e_1,e_2,\dots ,e_n$重复相同的过程, 从而倍数至多为$12n$.
再次地, 若Hessian稀疏且结构已知, 我们就能用图着色技术仅使用远少于$n$个的种子向量计算整个Hessian. 选取$q$的方式类似于有限差分中的选取方式. 而计算量的增加倍数至多为$12N_c({\nabla }^{2}f)$.

2.8 当前的限制

当前自动微分已经彰显了它在大型的和困难的优化问题中的能力. 但这样的工具却有可能在一些常用的优化框架和计算机计算时使人陷入困境.

1. 截断误差引发的导数危机. 例如, 若$f(x)$的计算依赖于求解偏微分方程(PDE), 则$f$的计算值会包含来自数值求解PDE时使用有限差分或有限元所引发的截断误差. 即我们有$\hat{f}(x)=f(x)+\tau (x)$, 其中$\hat{f}(\cdot )$为$f(\cdot )$的计算值, $\tau (\cdot )$为截断误差. 尽管$|\tau (x)|$通常来讲较小, 但它的导数${\tau }^{\prime }(x)$可就不一定了, 从而计算导数值$\hat{f}(\cdot )$中可能蕴含着巨大的误差(第1节中介绍的有限差分法也有相同的问题). 在计算机使用分段有理函数近似三角函数时, 也会出现类似的问题.
2. 代码分段导致计算失误. 另外一个问题的来源是, 将代码分段以提高特定区域上函数计算的速度或精确度. 一个病态的例子是$f(x)=x-1$. 若我们使用以下代码的分段计算该函数
   $$\mathrm{i}\mathrm{f}(x=1.0)\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}f=0.0\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}f=x-1.0,$$于是使用自动微分时我们就会有$f^{\prime }(1)=0$.

总的来说, 自动微分应当被视作一类愈发复杂的技术, 它们可以用来增强优化算法的性能, 从而使算法能够用于更广、更复杂的问题. 若提供敏度分析, 我们就能从计算的结果中挖掘更多的信息. 自动微分不是万能的. 我们不能因为有了自动微分就认为自己可以逃脱导数计算.

---

1. NP-completeness: https://en.wikipedia.org/wiki/NP-completeness ↩︎

2. G. N. NEWSAM AND J. D. RAMSDELL, Estimation of sparse Jacobian matrices, SIAM Journal on Algebraic and Discrete Methods, 4 (1983), pp. 404–418. ↩︎