机器学习中的自动微分

﹏ヽ暗。殇╰゛Y 2022-11-05 06:22 654阅读 0赞

### 目录 ###

*  摘要
 *  一 引言
 *  二 AD不是什么
 *   *  2.1 AD不是数值微分
     *  2.2 AD不是符号微分
 *  三. AD及其主要模式
 *   *  3.1 前向模式
     *   *  3.1.1 对偶数
     *  3.2 反向模式
     *  3.3 AD的起源和反向传播

--------------------

本文为文献Automatic Differentiation in Machine Learning: a Survey （2018）的原创译文，转载请注明出处。

# 摘要 #

导数，主要以梯度和Hessians的形式存在，在机器学习中随处可见。自动微分(AD)，也被称为算法微分或简单地称为“autodiff”，是一种与反向传播类似但比反向传播更通用的技术，用于有效和准确地计算以计算机程序表示的数值函数的导数。AD是一个小型但已建立起来的领域，其应用领域包括计算流体动力学、大气科学和工程设计优化。直到最近，机器学习和AD这两个领域在很大程度上还不了解彼此，而且在某些情况下，它们还各自独立地发现了彼此的研究结果。尽管它很重要，但通用AD并没有出现在机器学习工具箱中，随着它被称为“动态计算图”和“可微编程”，这种情况正在慢慢改变。我们调查了AD和机器学习的交集，涵盖了AD具有直接相关性的应用，并解决了主要应用技术。通过精确定义主要的微分技术及其相互关系，我们的目标是澄清术语“autodiff”、“automatic differentiation”和“符号微分”的用法，因为这些在机器学习中遇到的越来越多。

--------------------

# 一 引言 #

在计算机程序中微分方法可以大体被分为以下四类：

1.  手动推导计算微分并用程序实现；
2.  有限差值近似求出的数值微分；
3.  用表达式运算的符号微分（例如mathematica等）；
4.  自动微分，也被称为算法微分（本文重点）。

传统来说，机器学习中很多方法都用到大量的推导运算和依赖于对特定函数计算梯度和Hessians的传统算法，机器学习领域的研究者已经花费了很多精力试图将解析的导数手动推导出来并插入到像L-BFGS或随机梯度下降这样的标准的优化程序中。但是手动微分推导耗时且有误差，数值微分又会由于截断或约化等原因有较大误差，且对于梯度来说scales很差（？），使他不适用于机器学习这样需要对上百万参数进行求导的方法。符号微分可以解决手动和数值微分的这些劣势，但是却存在表达式过于复杂及表达式爆炸的问题。而且更重要的是，手动微分和符号微分都需要模型被定义为封闭形式的表达式，大大限制了算法控制性（？）和表达能力。  
我们得知第四种方法，自动微分（automatic diifferentiation）是非常强大的。自动微分其实是一种计算机程序的不正式的名字：这种计算机程序将变量的定义域替换为导数值和重新定义的沿微分链法则传播的导数运算符。除了其他广泛的运用外，通用的AD是在最近才被机器学习领域运用起来的。随着深度学习(LeCun等人2015年、Goodfellow等人2016年)在许多机器学习任务和现代工作流中处于最先进水平，基于快速原型和代码重用框架，如-Theano（Bastien等，2012），Torch（Collobert等，2011）和TensorFlow（Abadi等2016年)，情况正在慢慢改变，autograd2 (Maclaurin, 2016)、Chainer3 (Tokui等人，2015)和PyTorch4 (Paszke等人，2017)等项目正在引领通用AD进入主流。  
“自动”这个词可能会带来误导，使一些机器学习的实践者将自动微分或‘aotudiff’这样的标签用到任何不需要手动微分的地方，而不考虑其算法背后的数学原理。我们想强调AD是一个技术的术语，指的是一个特定的技术类别，这些技术通过在代码执行过程中数值导数的积累值而不是导数表达式来计算导数，从而生成数值导数计算结果。这允许导数值的精确度在机器精度只有很小的常数因数，并且允许导数值有理想化的渐近线效率（？）。与符号微分带来的编程上的麻烦相反，自动微分可以被用在改动最小的常规代码上，允许分支、循环和递归。由于这一普适性，AD被用于工业和学术的计算机模拟以及工程设计优化等领域(Forth和Evans, 2002;C’asanova等人，2002年)，计算流体动力学(Miller and c 'usdin, 2005: Thomas et al.， 2006;Bischof等人。2006)，物理建模(Ekstrom et al.， 2010)，最佳控制(Walther, 2007)，结构力学(Haase et al.， 2002)，大气科学(Carmichael and Sandu, 1997;Charpentier和Ghemires, 2000)，计算金融(Bischof等人，2002;Capriotti。2011)。  
在机器学习领域，与AD类似的一种专门的反向传播算法(backpropagation algorithm)一直是训练神经网络的主要手段，它在不同时期被独立研究人员重新发明，有着丰富多彩的历史(Griewank, 2012;Schmidhuber2015)。自从主要通过Rumelhart等人(1986)的工作而流行起来以来，它一直是研究和使用最多的训练算法之一。用最简单的术语来说，反向传播模型学习作为神经网络权值空间中的梯度下降，寻找目标函数的最小值。通过目标值在输出处的灵敏度反向传播(图1)得到所需的梯度，利用链法则计算目标对每个权重的偏导数。得到的算法本质上相当于将目标函数组成的网络评估函数在反向模式AD下进行变换，我们将看到，这实际上是对反向传播思想的推广。因此，对隐含于反向传播中的数学的适度理解为掌握AD技术提供了足够的背景。  
![图1][1]  
*图1：反向传播概览图。（a）训练输入 xi 被前向传播，得到相应的激活值 yi 。最终输出 y3 和目标输出 t 之间的差值 E 被计算出来。（b）差值的伴随值被反向传播，给出关于权重的梯度  ∇ w i E = ( ∂ E ∂ w 1 , . . . , ∂ E ∂ w 6 ) \\nabla\_\{w\_i\}E=(\\frac\{\\partial E\}\{\\partial w\_1\},...,\\frac\{\\partial E\}\{\\partial w\_6\}) ∇wiE=(∂w1∂E,...,∂w6∂E)，这也是随后被梯度下降法运用的计算步骤。关于输入值的导数 ∇ x i E \\nabla\_\{x\_i\} E ∇xiE也可以在同样的过程中计算出来。*

在本文中，我们从机器学习的角度回顾AD，包括它的起源，在机器学习中的应用，以及实现方法。在此过程中，我们还打算消除一些误解，我们认为这些误解阻碍了机器学习社区对AD的更广泛认识。在第二节中，我们首先说明AD与数字和符号的区别是如何不同的。第三节介绍了AD技术及其正、反累积方式。第4节讨论了衍生工具在机器学习中的作用，并审查了与AD相关的案例。第5节将介绍各种实现方法和通用的AD工具，接下来的第6节将讨论未来的方向。

--------------------

# 二 AD不是什么 #

没有合适的介绍，人们可能会认为AD是一种数值微分或符号微分。这个误解的产生是由于AD确实提供导数的数值（与表达式推导相反）并且它用到的是微分的符号规则（但与结果表达式不同的是它始终保留导数值），使他同时具有部分数值微分和部分符号微分的特点。

## 2.1 AD不是数值微分 ##

数值微分是用原始函数在特定样本点上的导数的有限差值近似（如图2）。在它的简洁的形式中，它基于导数有限的定义。例如对于一个多元函数 f : R n → R f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\} f:Rn→R，我们可以估计他的梯度 ∇ f = ( ∂ f ) ∂ x 1 , . . . , ∂ f ∂ x n ) \\nabla f = (\\frac\{\\partial f)\}\{\\partial x\_1\},...,\\frac\{\\partial f\}\{\\partial x\_n\}) ∇f=(∂x1∂f),...,∂xn∂f)用  
 ∂ f ( x ) ∂ x i ≈ f ( x + h e i ) − f ( x ) h ( 1 ) \\frac\{\\partial f(x)\}\{\\partial x\_i\} \\approx \\frac\{f(x+he\_i)-f(x)\}\{h\}\\quad \\quad (1) ∂xi∂f(x)≈hf(x\+hei)−f(x)(1)  
其中  e i e\_i ei 是第 i 个基矢，h > 0 是一个小的步长。这样是在实现起来不复杂上是有优势的，但是缺点在于对于 n 维导数计算 f 的 O(n) 是困难的，并且需要很认真的挑选步长。  
![图2][2]  
*图2：区分数学表达式和计算机代码的方法的范围，看一个截断的逻辑映射的例子(左上角)。符号微分(中右)给出了准确的结果，但需要封闭形式的输入，并存在表达式膨胀;数值微分(右下)由于舍入和截断误差而存在精度问题;自动微分(左下)与符号微分一样精确，只需要固定的常数因子和对控制流的支持。*

导数的数值估计从本质上就是不稳定的，除了应用于有限正定函数的复变方法。这是引入截断和舍入误差及步长的选取带来的。当  h → 0 h \\rightarrow0 h→0时截断误差趋于0。但是随着 h 减小，舍入误差逐渐增大并成为主要误差（如图3）。  
为了降低数值微分中的近似误差，人们发展了不同的技术，像中心差值近似  
 ∂ f ( x ) ∂ x i = f ( x + h e i ) − f ( x − h e i ) 2 h + O ( h 2 ) ( 2 ) \\frac\{\\partial f(x)\}\{\\partial x\_i\}=\\frac\{f(x+he\_i)-f(x-he\_i)\}\{2h\}+O(h^2) \\quad \\quad(2) ∂xi∂f(x)=2hf(x\+hei)−f(x−hei)\+O(h2)(2)  
其中一阶误差被取消，可以高效的将截断误差从一阶转移到 h 的二阶。对于一维的情况，计算（1）式的前向差值和（2）式的中心差值都只需要两次对 f 的计算。然而当维度增大，就将面临精确度和表现的矛盾，当计算一个  f : R n → R m f: \\mathbb\{R\}^n \\rightarrow \\mathbb\{R\}^m f:Rn→Rm的函数的雅可比矩阵就需要 2mn 计算。

![图3][3]  
*图3：（1）式中的前向误差和（2）式中的中心误差近似关于步长 h 的函数，截断逻辑图为  f ( x ) = 64 x ( 1 − x ) ( 1 − 2 x ) 2 ( 1 − 8 x + 8 x 2 ) 2 f(x)=64x(1-x)(1-2x)^2(1-8x+8x^2)^2 f(x)=64x(1−x)(1−2x)2(1−8x\+8x2)2。画出的误差是通过 E f o r w a r d ( h − x 0 ) = ∣ f ( x 0 + h ) − f ( x 0 ) h − d d x f ( x ) ∣ x 0 ∣ E\_\{forward\}(h- x\_0)=\\vert \\frac\{f(x\_0+h)-f(x\_0)\}\{h\}-\\frac\{d\}\{dx\}f(x)\\vert\_\{x\_0\}\\vert Eforward(h−x0)=∣hf(x0\+h)−f(x0)−dxdf(x)∣x0∣和 E c e n t e r ( h , x 0 ) = ∣ f ( x 0 + h ) − f ( x 0 − h ) 2 h − d d x f ( x ) ∣ x 0 ∣ E\_\{center\}(h,x\_0)=\\vert\\frac\{f(x\_0+h)-f(x\_0-h)\}\{2h\}-\\frac\{d\}\{dx\}f(x)\\vert\_\{x\_0\}\\vert Ecenter(h,x0)=∣2hf(x0\+h)−f(x0−h)−dxdf(x)∣x0∣ 在  x 0 = 0.2 x\_0=0.2 x0=0.2。*

还有一些其他的提高数值微分的方法，包括高阶有限差分、理查森极限外推法(Brezinski和Zaglia, 1991年)和使用加权和的微分求积法(Bert和Malik, 1996年)，都增加了计算的复杂性，并不能完全消除近似误差，仍然非常容易被浮点点截断。  
n维梯度的数值微分的O(n)复杂度是其在机器学习中应用的主要障碍，在最先进的深度学习模型中，n可能高达数百万或数十亿(Shazeer等人，2017)。相反，由于神经网络体系结构的错误复原能力得到了充分的证明，在深度学习设置中可以容忍近似错误(Gupta等人，2015)。

## 2.2 AD不是符号微分 ##

符号微分是一种针对表达式的自动操作，目的是得到微分表达式。运用微分法则例如  
 d d x ( f ( x ) + g ( x ) ) ⇝ d d x f ( x ) + d d x g ( x ) d d x ( f ( x ) g ( x ) ) ⇝ ( d d x f ( x ) ) g ( x ) + f ( x ) ( d d x g ( x ) ) ( 3 ) \\frac\{d\}\{dx\}(f(x)+g(x)) \\rightsquigarrow \\frac\{d\}\{dx\}f(x)+\\frac\{d\}\{dx\}g(x) \\\\ \\frac\{d\}\{dx\}(f(x)g(x)) \\rightsquigarrow (\\frac\{d\}\{dx\}f(x))g(x)+f(x)(\\frac\{d\}\{dx\}g(x)) \\quad (3) dxd(f(x)\+g(x))⇝dxdf(x)\+dxdg(x)dxd(f(x)g(x))⇝(dxdf(x))g(x)\+f(x)(dxdg(x))(3)

当公式被表示为数据结构时，用符号区分表达式树是一个完美的机械过程，甚至在微积分最开始的时候就被认为是机械自动化的对象(Leibniz, 1685)。这是在现代计算机代数系统，如Mathematica, Maxima，和Maple和机器学习框架，如Theano实现的。  
在优化中，符号微分可以为问题的定义域的结构提供有价值的见解，在某些情况下，可以产生极值的解析解(例如，求解 d d x f ( x ) = 0 \\frac\{d\}\{dx\}f(x) = 0 dxdf(x)=0)，从而完全消除导数计算的需要。另一方面，符号导数不能在运行时有效地计算导数值，因为它们可以得到比其导数所代表的表达式指数级增大了的表达式。  
考虑函数  h ( x ) = f ( x ) g ( x ) h(x)=f(x)g(x) h(x)=f(x)g(x)和（3）式中的乘法规则，由于 h 为乘积，  h ( x ) h(x) h(x)和 d d x h ( x ) \\frac\{d\}\{dx\}h(x) dxdh(x)都有共同的元素，  f ( x ) f(x) f(x)和 g ( x ) g(x) g(x)。在等式右侧， f ( x ) f(x) f(x)和 d d x f ( x ) \\frac\{d\}\{dx\}f(x) dxdf(x)单独出现。如果我们只是继续对 f ( x ) f(x) f(x)进行符号微分并将微分插入回合适的位置，我们就将嵌套进很多 f ( x ) f(x) f(x)和 d d x f ( x ) \\frac\{d\}\{dx\}f(x) dxdf(x)之间共同的重复项。因此，简单的符号微分会很容易导致指数级增多的符号表达式，从而导致运算时间很长。这个问题被称为表达式爆炸（如表1）。

*表1：逻辑图  l n + 1 = 4 l n ( 1 − l n ) ， l 1 = x l\_\{n+1\}=4l\_n(1-l\_n)，l\_1=x ln\+1=4ln(1−ln)，l1=x及 l n l\_n ln关于 x x x的导数的迭代，体现出表达式爆炸。*  
![表1][1 1]  
当我们考虑微分的精确数值估计并且符号形式不要太多，原则上通过仅在内存中存储中间的子表达式来大大简化计算是可行的。而且，为了更加高效，我们可以尽量多的交叉微分和简化步骤。这个交叉这两种步骤的思路形成了AD的基础，也为其提供了简单的模式：在初级操作层面应用符号微分并保留中间数值结果，与主函数的求值保持一致。这就是AD的前向积累模式，也是我们接下来要介绍的。

# 三. AD及其主要模式 #

AD可以被认为是对计算机程序的非标准解释，这种解释包括用各种导数的计算来增加标准计算。所有的数值计算最终都是由一组已知导数的有限初等运算组成的(Verma, 2000;Griewank和Walther, 2008)，通过链式法则结合各组成运算的导数，得出整体组成的导数。通常这些初等运算包括二元算术运算、一元符号转换和超越函数，如指数、对数和三角函数。  
在表2的左侧我们可以看到表达式 f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) f(x1,x2)=ln(x1)\+x1x2−sin(x2)的计算作为一个初等计算的计算轨迹被称为Wengert列表。我们使用Griewank 和 Walther三部分符号记录法，对于一个用中间变量 v i v\_i vi建立的函数 f : R n → R m f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\}^m f:Rn→Rm：

*  变量 v i − n = x i v\_\{i-n\}=x\_i vi−n=xi，i=1,…,n是输入变量；
 *  变量 v i , i = 1 , . . . , l v\_i,i=1,...,l vi,i=1,...,l为中间变量；
 *  变量 y m − i = v l − i , i = m − 1 , . . . , 0 y\_\{m-i\}=v\_\{l-i\},i=m-1,...,0 ym−i=vl−i,i=m−1,...,0为输出变量。

图4展示了给定初始操作的轨迹构成的计算图，有助于展示中间变量之间的依赖关系。  
计算的轨迹形成了AD算法的基础。这里一个重要的点是AD不仅能够对经典的封闭形式的表达式微分，还可以处理带有循环的算法，例如分支、循环、递归和调用过程等，这使它相比于符号微分具有了很大的优势。这要归功于对任意的数值代码最终都会成为一个带着确切的输入、中间或输出变量值的数值计算轨迹，这是我们用链法则计算微分时唯一要清楚的事情，不需要在意程序执行过程中的具体的控制流。另一个表达的方式是AD对任何操作都是会忽略的，包括控制流这种不会直接改变数值的操作。  
![图4][4]  
*图4：对于示例函数 f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) f(x1,x2)=ln(x1)\+x1x2−sin(x2)的计算图。*

*表2：AD前向模式示例。 y = f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) y=f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) y=f(x1,x2)=ln(x1)\+x1x2−sin(x2)在 ( x 1 , x 2 ) = ( 2 , 5 ) (x\_1,x\_2)=(2,5) (x1,x2)=(2,5)并设置 x ˙ 1 = 1 \\dot x\_1=1 x˙1=1来计算 ∂ y ∂ x 1 \\frac\{\\partial y\}\{\\partial x\_1\} ∂x1∂y。*  
![表2][2 1]

## 3.1 前向模式 ##

AD前向模式的积累是最简单的模式。考虑函数 f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) f(x1,x2)=ln(x1)\+x1x2−sin(x2)的计算轨迹如图4和表2。对于计算f关于 x 1 x\_1 x1的微分，我们从和每个中间变量 v i v\_i vi都有联系的微分开始  
 v ˙ i = ∂ v i ∂ x 1 \\dot v\_i = \\frac\{\\partial v\_i\}\{\\partial x\_1\} v˙i=∂x1∂vi  
对前向模式轨迹的每个基本操作都应用链法则，我们得到对应的微分轨迹，在表2中给出。在前后连接的步骤中计算原始值 v i v\_i vi和对应的导数值给出最终所要求的微分 v ˙ 5 = ∂ y ∂ x 1 \\dot v\_5=\\frac\{\\partial y\}\{\\partial x\_1\} v˙5=∂x1∂y。  
这自然地推广到计算有n个独立(输入)变量的函数 f : R n → R m f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\}^m f:Rn→Rm的雅可比矩阵;以及m个因变量 y j y\_j yj。在这种情况下，AD的每次向前传递都是通过只设置其中一个变量 x ˙ i = 1 \\dot x\_i=1 x˙i=1，并将其余的设置为0来初始化的(换句话说，设置 x ˙ i = e i \\dot x\_i = e\_i x˙i=ei，其中 e i e\_i ei为第i个单位向量)。使用指定的输入值 x = a x=a x=a运行代码，然后计算  
 y ˙ j = ∂ y j ∂ x i ∣ x = a , j = 1 , . . . , m , \\dot y\_j=\\frac\{\\partial y\_j\}\{\\partial x\_i\}\\vert\_\{x=a\},j=1,...,m, y˙j=∂xi∂yj∣x=a,j=1,...,m,  
给出雅可比矩阵的一列在a点的值。因此一个完整的雅可比矩阵能够通过n次计算得到。  
 J f = \[ ∂ y 1 ∂ x 1 . . . ∂ y 1 ∂ x n ⋮ ⋱ ⋮ ∂ y m ∂ x 1 . . . ∂ y m ∂ x n \] ∣ x = a \\mathbf\{J\_f\}= \\left\[ \\begin\{array\}\{ccc\} \\frac\{\\partial y\_1\}\{\\partial x\_1\} & ... & \\frac\{\\partial y\_1\}\{\\partial x\_n\}\\\\ \\vdots & \\ddots & \\vdots \\\\ \\frac\{\\partial y\_m\}\{\\partial x\_1\} & ... & \\frac\{\\partial y\_m\}\{\\partial x\_n\} \\end\{array\} \\right\] \\Bigg \\vert\_\{x=a\} Jf=⎣⎢⎡∂x1∂y1⋮∂x1∂ym...⋱...∂xn∂y1⋮∂xn∂ym⎦⎥⎤∣∣∣∣∣x=a  
另外，前向模式的AD提供了一个非常高效且不需要用到矩阵的方式，计算雅可比-向量乘积：  
 J f r = \[ ∂ y 1 ∂ x 1 . . . ∂ y 1 ∂ x n ⋮ ⋱ ⋮ ∂ y m ∂ x 1 . . . ∂ y m ∂ x n \] \[ r 1 ⋮ r n \] ( 4 ) \\mathbf\{J\_fr\}= \\left\[ \\begin\{array\}\{ccc\} \\frac\{\\partial y\_1\}\{\\partial x\_1\} & ... & \\frac\{\\partial y\_1\}\{\\partial x\_n\}\\\\ \\vdots & \\ddots & \\vdots \\\\ \\frac\{\\partial y\_m\}\{\\partial x\_1\} & ... & \\frac\{\\partial y\_m\}\{\\partial x\_n\} \\end\{array\} \\right\] \\left \[ \\begin\{array\}\{c\} r\_1\\\\ \\vdots \\\\ r\_n \\end\{array\} \\right\] \\quad (4) Jfr=⎣⎢⎡∂x1∂y1⋮∂x1∂ym...⋱...∂xn∂y1⋮∂xn∂ym⎦⎥⎤⎣⎢⎡r1⋮rn⎦⎥⎤(4)  
只需初始化 x ˙ = r \\dot x = r x˙=r。因此，我们可以只用一步前向过程来计算雅可比-向量乘积。对一个特殊的情况，当 f : R n → R f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\} f:Rn→R，我们可以得到关于给定向量 r r r的带方向的微分，用微分的线性组合来表示  
 ∇ f ⋅ r \\nabla f·r ∇f⋅r  
通过设置AD计算的初值 x ˙ = r \\dot x =r x˙=r。  
前向模式的AD对于 f : R → R m f:\\mathbb\{R\} \\rightarrow \\mathbb\{R\}^m f:R→Rm的函数是高效而直接的，因为所有的微分 d y i d x \\frac\{dy\_i\}\{dx\} dxdyi都能在一步前向传播中计算出来。相反的在另一种极限情况 f : R n → R f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\} f:Rn→R，前向模式AD就需要n次计算来得到梯度  
 ∇ f = ( ∂ y ∂ x 1 , . . . , ∂ y ∂ x n ， \\nabla f=(\\frac\{\\partial y\}\{\\partial x\_1\},...,\\frac\{\\partial y\}\{\\partial x\_n\}， ∇f=(∂x1∂y,...,∂xn∂y，  
这也对应着一个 1 × n 1 \\times n 1×n的雅可比矩阵在前向模式的n次运算中每次得到其中一列。  
总的来说，对于 f : R n → R m f:\\mathbb\{R\}^n \\rightarrow \\mathbb\{R\}^m f:Rn→Rm其中 n ≫ m n \\gg m n≫m的情况，我们会需要另一个方法来计算，我们将介绍反向模式的AD。

### 3.1.1 对偶数 ###

数学上看，前向模式的AD可以被看作用对偶数计算一个函数，可以被定义为阶段的泰勒展开的形式：  
 v + v ˙ ϵ v+\\dot v \\epsilon v\+v˙ϵ  
其中 v , v ˙ ∈ R v,\\dot v \\in \\mathbb\{R\} v,v˙∈R， ϵ \\epsilon ϵ是零幂数即 ϵ 2 = 0 , ϵ ≠ 0 \\epsilon^2=0,\\epsilon \\ne 0 ϵ2=0,ϵ=0。观察到，  
 ( v + v ˙ ϵ ) + ( u + u ˙ ϵ ) = ( v + u ) + ( v ˙ + u ˙ ) ϵ ( v + v ˙ ϵ ) ( u + u ˙ ϵ ) = ( v u ) + ( v u ˙ + v ˙ u ) ϵ (v+\\dot v\\epsilon)+(u+\\dot u\\epsilon)=(v+u)+(\\dot v+\\dot u)\\epsilon \\\\ (v+\\dot v\\epsilon)(u+\\dot u\\epsilon)=(vu)+(v\\dot u+\\dot v u)\\epsilon (v\+v˙ϵ)\+(u\+u˙ϵ)=(v\+u)\+(v˙\+u˙)ϵ(v\+v˙ϵ)(u\+u˙ϵ)=(vu)\+(vu˙\+v˙u)ϵ  
 ϵ \\epsilon ϵ的系数恰好是符号微分的反射（（3）式），我们可以通过设置  
 f ( v + v ˙ ϵ ) = f ( v ) + f ′ ( v ) v ˙ ϵ ( 5 ) f(v+\\dot v \\epsilon)=f(v)+f'(v)\\dot v \\epsilon \\quad \\quad (5) f(v\+v˙ϵ)=f(v)\+f′(v)v˙ϵ(5)  
来应用这一关系，并用对偶数作为数据结构来将导数值和原始值组合。链法则对于这一表达式也是适用的，（5）式的两个应用给出  
 f ( g ( v + v ˙ ϵ ) ) = f ( g ( v ) + g ′ ( v ) v ˙ ϵ ) = f ( g ( v ) ) + f ′ ( g ( v ) ) g ′ ( v ) v ˙ ϵ f(g(v+\\dot v\\epsilon))=f(g(v)+g'(v)\\dot v\\epsilon) \\\\ = f(g(v))+f'(g(v))g'(v)\\dot v\\epsilon f(g(v\+v˙ϵ))=f(g(v)\+g′(v)v˙ϵ)=f(g(v))\+f′(g(v))g′(v)v˙ϵ  
等式右侧 ϵ \\epsilon ϵ的系数恰好是  f f f和 g g g的组合的微分。这意味着我们关于不变的（5）式运用基本操作时所有的组成部分都要做，这反过来也意味着我们能够通过将任意非对偶数 v v v记作 v + 0 ϵ v+0\\epsilon v\+0ϵ来提取一个函数的微分并在这个不规范的方式取初值 ϵ \\epsilon ϵ的系数为1来计算函数值：  
 d f ( x ) d x ∣ x = v = ϵ 的 系 数 ( 对 偶 形 式 f ( v + 1 ϵ ) \\frac\{df(x)\}\{dx\}\\bigg \\vert\_\{x=v\}=\\epsilon的系数(对偶形式f(v+1\\epsilon) dxdf(x)∣∣∣∣x=v=ϵ的系数(对偶形式f(v\+1ϵ)  
这也扩展到任意的程序结构，因为对偶数作为数据类型，可以包含在任何数据结构中。只要一个对偶数保留在数据结构中，没有算术运算，它就会一直是一个对偶数;如果把它从数据结构中取出来，然后再进行操作，那么微分就会继续。  
在实践中，用选定的编程语言编写的函数f将被输入到AD工具中，然后AD工具将用相应的额外的代码对其进行扩充，以处理对偶操作，从而同时计算函数及其导数。这可以通过调用特定的库来实现，以源代码转换的形式，其中给定的源代码将被自动修改，或者通过操作符重载，使过程对用户透明。我们将在第5节讨论这些实现技术。

## 3.2 反向模式 ##

AD的反向积累模式对应于通常所说的反向传播算法，其中的微分传播时从给定的输出反向传播的。这是通过用一个伴随值来补充每个中间变量 v i v\_i vi：  
 v ˉ i = ∂ y j ∂ v i \\bar v\_i=\\frac\{\\partial y\_j\}\{\\partial v\_i\} vˉi=∂vi∂yj  
这代表了一个输出 y j y\_j yj关于 v i v\_i vi变化的敏感度。在反向传播下， y y y应为对应于误差E的标量。  
在反向模式的AD中，导数计算在两阶段过程的第二阶段。在第一阶段，原函数代码向前运行，填充中间变量 v i vi vi，并通过一个记录程序来记录计算图中的依赖关系。在第二阶段，通过传播伴随值 v ˉ i \\bar v\_i vˉi来计算导数，从输出到输入。  
回到这个例子 y = f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) y=f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) y=f(x1,x2)=ln(x1)\+x1x2−sin(x2)，在表3中，我们看到右边的伴随变量，对应于左边的每个原始初等运算。简单地说，我们感兴趣的是计算 v ˉ i = ∂ y ∂ v i \\bar v\_i=\\frac\{\\partial y\}\{\\partial v\_i\} vˉi=∂vi∂y关于每个变量 v i v\_i vi对于输出 y y y的贡献。以变量 v 0 v\_0 v0为例，从图4中可以看出，它影响 y y y的唯一方式是通过影响 v 2 v\_2 v2和 v 3 v\_3 v3，所以它对 y y y变化的贡献为  
 ∂ y ∂ v 0 = ∂ y ∂ v 2 ∂ v 2 ∂ v 0 + ∂ y ∂ v 3 ∂ v 3 ∂ v 0 o r v ˉ 0 = v ˉ 2 ∂ v 2 ∂ v 0 + v ˉ 3 ∂ v 3 ∂ v 0 \\frac\{\\partial y\}\{\\partial v\_0\}=\\frac\{\\partial y\}\{\\partial v\_2\}\\frac\{\\partial v\_2\}\{\\partial v\_0\}+\\frac\{\\partial y\}\{\\partial v\_3\}\\frac\{\\partial v\_3\}\{\\partial v\_0\} \\quad or \\quad \\bar v\_0=\\bar v\_2\\frac\{\\partial v\_2\}\{\\partial v\_0\}+\\bar v\_3\\frac\{\\partial v\_3\}\{\\partial v\_0\} ∂v0∂y=∂v2∂y∂v0∂v2\+∂v3∂y∂v0∂v3orvˉ0=vˉ2∂v0∂v2\+vˉ3∂v0∂v3  
在表3中，这一贡献被分为两步计算  
 v ˉ 0 = v ˉ 3 ∂ v 3 ∂ v 0 a n d v ˉ 0 = v ˉ 0 + v ˉ 2 ∂ v 2 ∂ v 0 \\bar v\_0=\\bar v\_3 \\frac\{\\partial v\_3\}\{\\partial v\_0\} \\quad and \\quad \\bar v\_0=\\bar v\_0+\\bar v\_2 \\frac\{\\partial v\_2\}\{\\partial v\_0\} vˉ0=vˉ3∂v0∂v3andvˉ0=vˉ0\+vˉ2∂v0∂v2  
将前向的轨迹中的线按他们原始的表达式排列了起来。  
在左手边向前传递之后，我们在右手边反向传递伴随项，从 v ˉ 5 = y ˉ = ∂ y ∂ y = 1 \\bar v\_5=\\bar y=\\frac\{\\partial y\}\{\\partial y\}=1 vˉ5=yˉ=∂y∂y=1 开始。最后我们在一次反向传递中就得到了导数  ∂ y ∂ x 1 = x 1 \\frac\{\\partial y\}\{\\partial x\_1\}=x\_1 ∂x1∂y=x1和 ∂ y ∂ x 2 = x ˉ 2 \\frac\{\\partial y\}\{\\partial x\_2\}=\\bar x\_2 ∂x2∂y=xˉ2。  
与正向积累模式的直接性相比，反向模式的AD起初可能显得有些“神秘”(Dennis and Schnabel, 1996)。Griewank （2008)认为，这一部分是因为人们普遍熟悉链式法则，将其作为了一种向前传播导数的机械程序。  
反向模式的一个重要优点是，对于具有大量输入的函数，反向模式的求值(就操作计数而言)成本明显低于前向模式。在 f : R n → R f: \\mathbb\{R\}^n→\\mathbb\{R\} f:Rn→R的极端情况下，只要一步反向模式就足以计算全部的梯度 ∇ f = ( ∂ y ∂ x 1 , . . . , ∂ y ∂ x n ) \\nabla f=(\\frac\{\\partial y\}\{\\partial x\_1\},...,\\frac\{\\partial y\}\{\\partial x\_n\}) ∇f=(∂x1∂y,...,∂xn∂y)，而填充相同的数据则需要n次正向模式。由于机器学习实践主要涉及标量值目标相对于大量参数的梯度，这建立了反向模式，而不是以反向传播算法的形式作为支柱的正向模式。

*表3：反向模式AD的例子。计算 y = f ( x 1 , x 2 ) = l n ( x 1 ) + x 1 x 2 − s i n ( x 2 ) y=f(x\_1,x\_2)=ln(x\_1)+x\_1x\_2-sin(x\_2) y=f(x1,x2)=ln(x1)\+x1x2−sin(x2)在 ( x 1 , x 2 ) = ( 2 , 5 ) (x\_1,x\_2)=(2,5) (x1,x2)=(2,5)处的值。在左侧前向模式计算初始值后，右侧伴随操作被反向计算。注意 ∂ y ∂ x 1 \\frac\{\\partial y\}\{\\partial x\_1\} ∂x1∂y和 ∂ y ∂ x 2 \\frac\{\\partial y\}\{\\partial x\_2\} ∂x2∂y在反向步骤中是同时被计算的，从伴随值 v ˉ 5 = y ˉ = ∂ y ∂ y = 1 \\bar v\_5=\\bar y=\\frac\{\\partial y\}\{\\partial y\}=1 vˉ5=yˉ=∂y∂y=1。*  
![表3][3 1]  
一般来说,对于一个函数 f : R n → R m f: \\mathbb\{R\}^n \\rightarrow \\mathbb\{R\}^m f:Rn→Rm，如果我们将用来评估最初函数的操作数记为 o p s ( f ) ops(f) ops(f)，通过前向模式计算 m × n m \\times n m×n雅可比矩阵的时间记为 n c o p s ( f ) n c ops(f) ncops(f),而相同的计算反向模式可以在 m c o p s ( f ) m c ops(f) mcops(f)的时间内完成,其中 c是一个常数， c < 6 c <6 c<6 ,  c   \[ 2 , 3 \] c ~\[2, 3\] c \[2,3\]（Griewank沃尔特,2008)。也就是说，反向模式AD在 m ≪ n m \\ll n m≪n时表现更好。

与正向模式下雅可比向量积的无矩阵计算(式4)类似，反向模式可以用于转置的Jacobian-vector乘积计算：  
 J f T r = \[ ∂ y 1 ∂ x 1 . . . ∂ y m ∂ x 1 ⋮ ⋱ ⋮ ∂ y 1 ∂ x n . . . ∂ y m ∂ x n \] \[ r 1 ⋮ r m \] \\mathbf\{J^T\_fr\}= \\left\[ \\begin\{array\}\{ccc\} \\frac\{\\partial y\_1\}\{\\partial x\_1\} & ... & \\frac\{\\partial y\_m\}\{\\partial x\_1\}\\\\ \\vdots & \\ddots & \\vdots \\\\ \\frac\{\\partial y\_1\}\{\\partial x\_n\} & ... & \\frac\{\\partial y\_m\}\{\\partial x\_n\} \\end\{array\} \\right\] \\left \[ \\begin\{array\}\{c\} r\_1\\\\ \\vdots \\\\ r\_m \\end\{array\} \\right\] JfTr=⎣⎢⎡∂x1∂y1⋮∂xn∂y1...⋱...∂x1∂ym⋮∂xn∂ym⎦⎥⎤⎣⎢⎡r1⋮rm⎦⎥⎤

通过初始化反向模式 y ˉ = r \\bar y = r yˉ=r  
然而，反向模式AD的优点随着存储需求成本的增加(在最坏的情况下)，存储需求的增长会与所评估函数中的操作数量成正比。通过使用检查点策略和数据流分析等先进方法来提高实现中的存储需求是一个活跃的研究领域(Dauvergne和Hascoêt, 2006;Siskind and Pearlmutter, 2017)。

## 3.3 AD的起源和反向传播 ##

AD背后的理念可以追溯到20世纪50年代(诺兰，1953年;Beda等人，1959年)。正向模式AD作为求偏导数的一般方法，本质上是由Wengert(1964)发现的。随后是一段活动相对较低的时期，直到1980年代主要通过Griewank(1989年)的工作重新引起人们对该领域的兴趣，现代程序设计语言的改进和有效的反向模式广告的可行性也提供了支持。反向模式广告和反向传播有着交织的历史。在连续时间形式主义中，逆向模式的本质是庞特里亚金最大原理(Rozonoer, 1959;boltyanski等人，1960年)。这种方法被控制理论团体所理解(Bryson and Denham, 1962;(Bryson and Ho, 1969)，并通过Werbos(1974)将离散时间变量按依赖关系拓扑排序，以更正式的形式进行转换。在Werbos之前，Linnainmaa(1970, 1976)的著作经常被引用为第一个描述逆向模式的著作。Speelpenning(1980)随后引入了逆向模式AD，从这个意义上说，他给出了第一个实际上是自动的实现，接受用通用编程语言编写的计算过程规范，并自动执行逆向模式转换。顺便提一句，Hecht-Nielsen(1989)引用了Bryson and Ho(1969)和Werbos(1974)的工作作为已知的最早的反向传播实例。在机器学习社区中，该方法已经被重新发明了好几次，例如Parker(1985)，直到它最终被Rumelhart等人(1986)和并行分布式处理(PDP)小组声名鹊起。PDP小组是在他们自己的发现之后才知道帕克的工作的;同样，Werbos的工作直到被Parker发现才得到赞赏(Hecht-Nielsen, 1989)。这给我们讲了一个有趣的故事:两个高度相互关联的研究团体却在这个奠基时期设法保持了分离。

--------------------

[1]: /images/20221023/ff0f65ff94f94a20bc30b9292c24b0c6.png
[2]: /images/20221023/8918e70a8fd5415fb2ae089c350ff444.png
[3]: /images/20221023/7049f128024743df9edf2d4e5175dc97.png
[1 1]: /images/20221023/2132b29d936843ff96903a7211d0fceb.png
[4]: /images/20221023/d34abe2659654cf884698b7b51a56e2a.png
[2 1]: /images/20221023/3f49242bde8e4f2f8bf873baa6cbac2c.png
[3 1]: /images/20221023/f52de3dec66e4802b603c55eccb29daf.png