(转)自动微分(Automatic Differentiation)简介——tensorflow核心原理

现代深度学习系统中(好比MXNet, TensorFlow等)都用到了一种技术——自动微分。在此以前,机器学习社区中不多发挥这个利器,通常都是用Backpropagation进行梯度求解,而后进行SGD等进行优化更新。手动实现过backprop算法的同窗应该能够体会到其中的复杂性和易错性,一个好的框架应该能够很好地将这部分难点隐藏于用户视角,而自动微分技术刚好能够优雅解决这个问题。接下来咱们将一块儿学习这个优雅的技术:-)。本文主要来源于陈天奇在华盛顿任教的课程CSE599G1: Deep Learning System和《Automatic differentiation in machine learning: a survey》。css

什么是自动微分

微分求解大体能够分为4种方式:node

  • 手动求解法(Manual Differentiation)
  • 数值微分法(Numerical Differentiation)
  • 符号微分法(Symbolic Differentiation)
  • 自动微分法(Automatic Differentiation)

为了讲明白什么是自动微分,咱们有必要了解其余方法,作到有区分有对比,从而更加深刻理解自动微分技术。git

手动求解法

手动求解其实就对应咱们传统的backprop算法,咱们求解出梯度公式,而后编写代码,代入实际数值,得出真实的梯度。在这样的方式下,每一次咱们修改算法模型,都要修改对应的梯度求解算法,所以没有很好的办法解脱用户手动编写梯度求解的代码,这也是为何咱们须要自动微分技术的缘由。github

数值微分法

数值微分法是根据导数的原始定义: 算法

f′(x)=limh→0f(x+h)−f(x)hf′(x)=limh→0f(x+h)−f(x)hexpress


那么只要hh取很小的数值,好比0.0001,那么咱们能够很方便求解导数,而且能够对用户隐藏求解过程,用户只要给出目标函数和要求解的梯度的变量,程序能够自动给出相应的梯度,这也是某种意义上的“自动微分”:-)。不幸的是,数值微分法计算量太大,求解速度是这四种方法中最慢的,更加雪上加霜的是,它引发的roundoff errortruncation error使其更加不具有实际应用场景,为了弥补缺点,便有以下center difference approximation: 编程

f′(x)=limh→0f(x+h)−f(x−h)2hf′(x)=limh→0f(x+h)−f(x−h)2h数组


惋惜并不能彻底消除truncation error,只是将偏差减少。虽然数值微分法有如上缺点,可是因为它实在是太简单实现了,因而不少时候,咱们利用它来检验其余算法的正确性,好比在实现backprop的时候,咱们用的”gradient check”就是利用数值微分法。缓存

 

符号微分法

符号微分是代替咱们第一种手动求解法的过程,利用代数软件,实现微分的一些公式好比: 网络

ddx(f(x)+g(x))=ddxf(x)+ddxg(x)ddxf(x)g(x)=(ddxf(x))g(x)+f(x)(ddxg(x))ddxf(x)g(x)=f′(x)g(x)−f(x)g′(x)g(x)2ddx(f(x)+g(x))=ddxf(x)+ddxg(x)ddxf(x)g(x)=(ddxf(x))g(x)+f(x)(ddxg(x))ddxf(x)g(x)=f′(x)g(x)−f(x)g′(x)g(x)2


而后对用户提供的具备closed form的数学表达式进行“自动微分”求解,什么是具备closed form的呢?也就是必须能写成完整数学表达式的,不能有编程语言中的循环结构,条件结构等。所以若是能将问题转化为一个纯数学符号问题,咱们能利用现有的代数软件进行符号微分求解,这种程度意义上的“自动微分”其实已经很完美了。然而缺点咱们刚刚也说起过了,就是必需要有closed form的数学表达式,另外一个有名的缺点是“表达式膨胀”(expression swell)问题,若是不加当心就会使得问题符号微分求解的表达式急速“膨胀”,致使最终求解速度变慢,对于这个问题请看以下图: 
这里写图片描述
稍不注意,符号微分求解就会如上中间列所示,表达式急剧膨胀,致使问题求解也随着变慢。

 

自动微分法

终于轮到咱们的主角登场,自动微分的存在依赖于它识破以下事实:

全部数值计算归根结底是一系列有限的可微算子的组合

自动微分法是一种介于符号微分和数值微分的方法:数值微分强调一开始直接代入数值近似求解;符号微分强调直接对代数进行求解,最后才代入问题数值;自动微分将符号微分法应用于最基本的算子,好比常数,幂函数,指数函数,对数函数,三角函数等,而后代入数值,保留中间结果,最后再应用于整个函数。所以它应用至关灵活,能够作到彻底向用户隐藏微分求解过程,因为它只对基本函数或常数运用符号微分法则,因此它能够灵活结合编程语言的循环结构,条件结构等,使用自动微分和不使用自动微分对代码整体改动很是小,而且因为它的计算实际是一种图计算,能够对其作不少优化,这也是为何该方法在现代深度学习系统中得以普遍应用。

自动微分Forward Mode

考察以下函数: 

f(x1,x2)=ln(x1)+x1x2−sin(x2)f(x1,x2)=ln(x1)+x1x2−sin(x2)


咱们能够将其转化为以下计算图: 
这里写图片描述
转化成如上DAG(有向无环图)结构以后,咱们能够很容易分步计算函数的值,并求取它每一步的导数值: 
这里写图片描述
上表中左半部分是从左往右每一个图节点的求值结果,右半部分是每一个节点对于x1x1的求导结果,好比v1˙=dvdx1v1˙=dvdx1,注意到每一步的求导都利用到上一步的求导结果,这样不至于重复计算,所以也不会产生像符号微分法的”expression swell”问题。 
自动微分的forward mode很是符合咱们高数里面学习的求导过程,只要您对求导法则还有印象,理解forward mode自不在话下。若是函数输入输出为: 

R→RmR→Rm


那么利用forward mode只需计算一次如上表右边过程便可,很是高效。对于输入输出映射为以下的: 

Rn→RmRn→Rm


这样一个有nn个输入的函数,求解函数梯度须要nn遍如上计算过程。然而实际算法模型中,好比神经网络,一般输入输出是极其不成比例的,也就是: 

n>>mn>>m


那么利用forward mode进行自动微分就过低效了,所以便有下面要介绍的reverse mode。

 

自动微分Reverse Mode

若是您理解神经网络的backprop算法,那么恭喜你,自动微分的backward mode其实就是一种通用的backprop算法,也就是backprop是reverse mode自动微分的一种特殊形式。从名字能够看出,reverse mode和forward mode是一对相反过程,reverse mode从最终结果开始求导,利用最终输出对每个节点进行求导,其过程以下红色箭头所示: 
这里写图片描述 
其具体计算过程以下表所示: 
这里写图片描述
上表左边和以前的forward mode一致,用于求解函数值,右边则是reverse mode的计算过程,注意必须从下网上看,也就是一开始先计算输出yy对于节点v5v5的导数,用v¯¯¯5v¯5表示dydv5dydv5,这样的记号能够强调咱们对当前计算结果进行缓存,以便用于后续计算,而没必要重复计算。由链式法则咱们能够计算输出对于每一个节点的导数。 
好比对于节点v3v3: 

dydv3=dydv5dv5dv3dydv3=dydv5dv5dv3


用另外一种记法变获得: 

dydv3=v5¯¯¯¯¯dv5dv3dydv3=v5¯dv5dv3


好比对于节点v0v0: 

dydv0=dydv2dv2dv0+dydv3dv3dv0dydv0=dydv2dv2dv0+dydv3dv3dv0


若是用另外一种记法,即可得出: 

dydv0=v¯¯¯2dv2dv0+v¯¯¯3dv3dv0dydv0=v¯2dv2dv0+v¯3dv3dv0


和backprop算法同样,咱们必须记住前向时当前节点发出的边,而后在后向传播的时候,能够搜集全部受到当前节点影响节点。 
如上的计算过程,对于像神经网络这种模型,一般输入是上万到上百万维,而输出损失函数是1维的模型,只须要一遍reverse mode的计算过程,即可以求出输出对于各个输入的导数,从而轻松求取梯度用于后续优化更新。

 

自动微分的实现

这里主要讲解reverse mode的实现方式,forward mode的实现基本和reverse mode一致,可是因为机器学习算法中大部分用reverse mode才能够高效求解,因此它是咱们理解的重心。代码设计轮廓来源于CSE599G1的做业,经过分析完成做业,能够展现自动微分的简洁性和灵活可用性。 
首先自动微分会将问题转化成一种有向无环图,所以咱们必须构造基本的图部件,包括节点和边。能够先看看节点是如何实现的: 
这里写图片描述
首先节点能够分为三种:

  • 常数节点
  • 变量节点
  • 带操做算子节点

所以Node类中定义了op成员用于存储节点的操做算子,const_attr表明节点的常数值,name是节点的标识,主要用于调试。 
对于边的实现则简单的多,每一个节点只要知道自己的输入节点便可,所以用inputs来描述节点的关系。 
有了如上的定义,利用操做符重载,咱们能够很简单构造一个计算图,举一个简单的例子: 

f(x1,x2)=x1x2+x2f(x1,x2)=x1x2+x2


对于如上函数,只要重载加法和乘法操做符,咱们能够轻松获得以下计算图: 
这里写图片描述

 

操做算子是自动微分最重要的组成部分,接下来咱们重点介绍,先上代码: 
这里写图片描述
从定义能够看出,全部实际计算都落在各个操做算子中,上面代码应该抽象一些,咱们来举一个乘法算子的例子加以说明: 
这里写图片描述
咱们重点讲解一下gradient方法,它接收两个参数,一个是node,也就是当前要计算的节点,而output_grad则是后面节点传来的,咱们来看看它究竟是啥玩意,对于以下例子: 

y=f(x1∗x2)y=f(x1∗x2)


那么要求yy关于x1x1的导数,那么根据链式法则可得: 

∂y∂x1=∂y∂f∂f∂x1=∂y∂x1x2∂x1x2∂x1=output_grad∗x2∂y∂x1=∂y∂f∂f∂x1=∂y∂x1x2∂x1x2∂x1=output_grad∗x2


则output_grad就是上面的∂y∂f∂y∂f,计算yy对于x2x2相似。所以在程序中咱们会返回以下:

 

return [node.inputs[1] * output_grad, node.inputs[0] * output_grad]
  • 1

再来介绍一个特殊的op——PlaceHolderOp,它的做用就如同名字,起到占位符的做用,也就是自动微分中的变量,它不会参与实际计算,只等待用户给他提供实际值,所以他的实现以下: 
这里写图片描述

了解了节点和操做算子的定义,接下来咱们考虑如何协调执行运算。首先是如何计算函数值,对于一幅计算图,因为节点与节点之间的计算有必定的依赖关系,好比必须先计算node1以后才能够计算node2,那么如何能正确处理好计算关系呢?一个简单的方式是对图节点进行拓扑排序,这样能够保证须要先计算的节点先获得计算。这部分代码由Executor掌控: 
这里写图片描述

Executor是实际计算图的引擎,用户提供须要计算的图和实际输入,Executor计算相应的值和梯度。

如何从计算图中计算函数的值,上面咱们已经介绍了,接下来是如何自动计算梯度。reverse mode的自动微分,要求从输出到输入节点,按照前后依赖关系,对各个节点求取输出对于当前节点的梯度,那么和咱们上面介绍的恰好相反,为了获得正确计算节点顺序,咱们能够将图节点的拓扑排序倒序便可。代码也很简单,以下所示: 
这里写图片描述
这里先介绍一个新的算子——oneslike_op。他是一个和numpy自带的oneslike函数同样的算子,做用是构造reverse梯度图的起点,由于最终输出关于自己的梯度就是一个和输出shape同样的全1数组,引入oneslike_op可使得真实计算得以延后,所以gradients方法最终返回的不是真实的梯度,而是梯度计算图,而后能够复用Executor,计算实际的梯度值。 
紧接着是根据输出节点,得到倒序的拓扑排序序列,而后遍历序列,构造实际的梯度计算图。咱们重点来介绍node_to_output_grad和node_to_output_grads_list这两个字典的意义。 
先关注node_to_output_grads_list,他key是节点,value是一个梯度列表,表明什么含义呢?先看以下部分计算图: 
这里写图片描述 
此时咱们要计算输出yy关于节点n1n1的导数,那么咱们观察到他的发射边链接的节点有n3,n4n3,n4,而对应n3,n4n3,n4节点调用相应op的gradient方法,会返回输出yy关于各个输入节点的导数。此时为了准确计算输出yy关于节点n1n1的导数,咱们须要将其发射边关联节点的计算梯度搜集起来,好比上面的例子,咱们须要搜集: 

node_to_output_grads_list={n1:[∂y∂n3∂n3∂n1,∂y∂n4∂n4∂n1]}node_to_output_grads_list={n1:[∂y∂n3∂n3∂n1,∂y∂n4∂n4∂n1]}


一旦搜集好对应输出边节点关于当前节点导数,那么当前节点的导数即可以由链式法则计算得出,也就是: 

∂y∂n1=∂y∂n3∂n3∂n1+∂y∂n4∂n4∂n1∂y∂n1=∂y∂n3∂n3∂n1+∂y∂n4∂n4∂n1


所以node_to_output_grad字典存储的就是节点对应的输出关于节点的导数。通过gradients函数执行后,会返回须要求取输出关于某节点的梯度计算图: 
这里写图片描述 
而对于Executor而言,它并不知道此时的图是否被反转,它只关注用户实际输入,还有计算相应的值而已。

 

自动梯度的应用

有了上面的大篇幅介绍,咱们其实已经实现了一个简单的自动微分引擎了,接下来看如何使用: 
这里写图片描述
使用至关简单,咱们像编写普通程序同样,对变量进行各类操做,只要提供要求导数的变量,还有提供实际输入,引擎能够正确给出相应的梯度值。 
下面给出一个根据自动微分训练Logistic Regression的例子:

  1.  
    import autodiff as ad
  2.  
    import numpy as np
  3.  
     
  4.  
     
  5.  
    def logistic_prob(_w):
  6.  
    def wrapper(_x):
  7.  
    return 1 / (1 + np.exp(-np.sum(_x * _w)))
  8.  
    return wrapper
  9.  
     
  10.  
     
  11.  
    def test_accuracy(_w, _X, _Y):
  12.  
    prob = logistic_prob(_w)
  13.  
    correct = 0
  14.  
    total = len(_Y)
  15.  
    for i in range(len(_Y)):
  16.  
    x = _X[i]
  17.  
    y = _Y[i]
  18.  
    p = prob(x)
  19.  
    if p >= 0.5 and y == 1.0:
  20.  
    correct += 1
  21.  
    elif p < 0.5 and y == 0.0:
  22.  
    correct += 1
  23.  
    print( "总数:%d, 预测正确:%d" % (total, correct))
  24.  
     
  25.  
     
  26.  
    def plot(N, X_val, Y_val, w, with_boundary=False):
  27.  
    import matplotlib.pyplot as plt
  28.  
    for i in range(N):
  29.  
    __x = X_val[i]
  30.  
    if Y_val[i] == 1:
  31.  
    plt.plot(__x[ 1], __x[2], marker='x')
  32.  
    else:
  33.  
    plt.plot(__x[ 1], __x[2], marker='o')
  34.  
    if with_boundary:
  35.  
    min_x1 = min(X_val[:, 1])
  36.  
    max_x1 = max(X_val[:, 1])
  37.  
    min_x2 = float(-w[ 0] - w[1] * min_x1) / w[2]
  38.  
    max_x2 = float(-w[ 0] - w[1] * max_x1) / w[2]
  39.  
    plt.plot([min_x1, max_x1], [min_x2, max_x2], '-r')
  40.  
     
  41.  
    plt.show()
  42.  
     
  43.  
     
  44.  
    def gen_2d_data(n):
  45.  
    x_data = np.random.random([n, 2])
  46.  
    y_data = np.ones(n)
  47.  
    for i in range(n):
  48.  
    d = x_data[i]
  49.  
    if d[0] + d[1] < 1:
  50.  
    y_data[i] = 0
  51.  
    x_data_with_bias = np.ones([n, 3])
  52.  
    x_data_with_bias[:, 1:] = x_data
  53.  
    return x_data_with_bias, y_data
  54.  
     
  55.  
     
  56.  
    def auto_diff_lr():
  57.  
    x = ad.Variable(name= 'x')
  58.  
    w = ad.Variable(name= 'w')
  59.  
    y = ad.Variable(name= 'y')
  60.  
     
  61.  
    # 注意,如下实现某些状况会有很大的数值偏差,
  62.  
    # 因此通常真实系统实现会提供高阶算子,从而减小数值偏差
  63.  
     
  64.  
    h = 1 / (1 + ad.exp(-ad.reduce_sum(w * x)))
  65.  
    L = y * ad.log(h) + ( 1 - y) * ad.log(1 - h)
  66.  
    w_grad, = ad.gradients(L, [w])
  67.  
    executor = ad.Executor([L, w_grad])
  68.  
     
  69.  
    N = 100
  70.  
    X_val, Y_val = gen_2d_data(N)
  71.  
    w_val = np.ones( 3)
  72.  
     
  73.  
    plot(N, X_val, Y_val, w_val)
  74.  
    executor = ad.Executor([L, w_grad])
  75.  
    test_accuracy(w_val, X_val, Y_val)
  76.  
    alpha = 0.01
  77.  
    max_iters = 300
  78.  
    for iteration in range(max_iters):
  79.  
    acc_L_val = 0
  80.  
    for i in range(N):
  81.  
    x_val = X_val[i]
  82.  
    y_val = np.array(Y_val[i])
  83.  
    L_val, w_grad_val = executor.run(feed_dict={w: w_val, x: x_val, y: y_val})
  84.  
    w_val += alpha * w_grad_val
  85.  
    acc_L_val += L_val
  86.  
    print( "iter = %d, likelihood = %s, w = %s" % (iteration, acc_L_val, w_val))
  87.  
    test_accuracy(w_val, X_val, Y_val)
  88.  
    plot(N, X_val, Y_val, w_val, True)
  89.  
     
  90.  
     
  91.  
    if __name__ == '__main__':
  92.  
    auto_diff_lr()
相关文章
相关标签/搜索