前言

这篇文章先介绍一下自动微分吧，这一个月来对我来说真的是多事之秋，能遇上的各种事情几乎都遇上了，但愿生活节奏和工作方面能够慢慢平稳下来。

第一部分的传送门:机器学习系统或者SysML&DL笔记(一)，第二部分是要介绍一些网络结构，但是没有完善，所以先把第三部分先发出来。

接下来我们就正式开始进入深度学习库的内部，来探索其具体的实现过程了。这节课的主要话题是Backpropagation and Automatic Differentiation，对于深度学习从业者，这是一个必须要掌握的技术，我们不光要理解其原理，更多的是将代码手撸出来。

以下是本系列文章的笔记来源：

• CSE 599W: Systems for ML
• AI-Sys Spring 2019

同时本文部分内容也参考深度学习花书第六章《深度前馈网络》。

微分方法

我们在实际训练中，模型训练的过程主要分这几步：

• 输入图像
• 图像通过神经网络层提取到了特征信息
• 将网络层最后输出的信息输入predictor中产生最终的预测结果，一般神经网络的最后一层为输出层，这个例子中的输出层就是sigmoid层
• 根据预测结果得到损失值，进行训练(不同任务采取的损失函数和优化器不同)

正如上图所示，使用的predictor为Sigmoid function：
$S(x)=\frac{1}{1+e^{-x}}=\frac{e^{x}}{e^{x}+1}$

这个函数可以将输入的数据转化0-1之间，充当激活神经元的作用。那么我们既然得到了神经网络的输出，接下来就是将输出与标签信息一同输入损失函数($L(w)$)去计算出损失值，随后通过优化器(SGD)对权重数据($w$)进行更新。

上述这个过程中涉及到了两个子过程：

• 前向传播(前馈网络可以被视为一种高效的非线性函数近似器)
• 反向求导与梯度更新(仅仅使得损失函数达到一个非常小的值，但并不能保证全局收敛)

线性模型，如逻辑回归和线性回归，是非常吸引人的 ，因为无论是通过闭解形式还是还是用凸优化，他都能够高效且可靠地拟合。线性模型也有明显的缺陷，那就是该模型的能力被局限在线性函数里，所以它无法理解任何两个输入变量见的相互作用

好了稍微有点啰嗦，那么这个过程中，反向求导的计算部分是怎样的？我们如何求导去更新梯度的值呢？这个时候就涉及到了自动微分的概念，说起自动微分，我们先列举下目前存在的一些微分方法：

• 手动微分
• 数值微分
• 符号微分
• 自动微分

手动微分

手动微分，顾名思义，手动微分法就是需要我们手动编写出代价函数、激活函数的求导代码，硬编码这些函数的求导方法，如果这些函数后面有调整该函数的求导方法又要重新实现，可以说是又麻烦又容易出错。

那为什么还需要手动微分，是因为有些我们自定义的算子可能比较复杂，涉及到了新的运算，目前深度学习库中实现的算子还没有提供。或者说这个算子比较特殊，虽然可以用很多算子组合起来，但是速度方面没有直接重新写一个一体的高效率，这样下来就需要重新手动写前向和反向过程了。

例如在Pytorch中我们可以使用C++拓展Cuda来编写前向和反向求导功能函数接入Pytorch中当成拓展使用。

数值微分

数值微分就是利用导数的定义来实现的，我们只需要一个F函数(目标函数)和一个需要求解的输入x即可计算出相关的梯度。当h(一般取值1e-5或者1e-6)取值很小的时候，下方的等式成立：

但是这种计算方式有两个问题(自行维基百科)：

• Roundoff Error
• Truncation Error

下面的一种使用Center形式计算导数的方法可以缓解Truncation错误但是依然会有可能Round的错误。

总结一句，不具备现实的使用场景，一般用于检验其他方式求解出来导数的正确性，例如在CS231n中notebook中经常遇到的：

# As we did for the SVM, use numeric gradient checking as a debugging tool.
# The numeric gradient should be close to the analytic gradient.
from cs231n.gradient_check import grad_check_sparse
f = lambda w: softmax_loss_naive(w, X_dev, y_dev, 0.0)[0]
grad_numerical = grad_check_sparse(f, W, grad, 10)

其中grad_check_sparse的具体实现为：

def grad_check_sparse(f, x, analytic_grad, num_checks=10, h=1e-5):
    """
    sample a few random elements and only return numerical
    in this dimensions.
    """

    for i in range(num_checks):
        ix = tuple([randrange(m) for m in x.shape])

        oldval = x[ix]
        x[ix] = oldval + h # increment by h
        fxph = f(x) # evaluate f(x + h)
        x[ix] = oldval - h # increment by h
        fxmh = f(x) # evaluate f(x - h)
        x[ix] = oldval # reset

        grad_numerical = (fxph - fxmh) / (2 * h)
        grad_analytic = analytic_grad[ix]
        rel_error = (abs(grad_numerical - grad_analytic) /
                    (abs(grad_numerical) + abs(grad_analytic)))
        print('numerical: %f analytic: %f, relative error: %e'
              %(grad_numerical, grad_analytic, rel_error))

符号微分

符号微分作为一种比较常用的微分法，在Matlab、Octave软件中我们经常能够遇到。

符号微分，顾名思义就是利用符号上面的算式来对微分符号表达式进行简化，然后进行求解。此时我们输入的是一个用符号表达式组成的树(或者叫做计算图)，然后通过微分规则对这些表达式进行化简。使用到的微分规则有sum rule、product rule、chain rule等等，也就是利用计算机去帮你微分，也因此限定了规则，必须是closed-form，不能有循环和条件结构等等。

这种方法一旦遇到非常复杂的公式很容易导致“表达式膨胀”，就如上图中间的公式所示，这种情况下计算速度会大大降低。

总结下，虽然这样来说，计算梯度的解析表达式是很直观的，但是数值化地求解这样的表示是在计算上的代价可能会很大。

BP(BackproPagation) 反向传播

说道自动微分，必须提一下反向传播。

反向传播是我们听过最多的一种反向求导的方法，在CS231n中有较为详细的讲解，所以这里不进行详细的描述了，仅仅是说一下反向求导的缺点：

• 我们需要保存前向传播中的中间输出变量以在反向求导的时候使用它
• 缺少灵活性，例如我们会计算梯度的梯度(gradient of gradient)。

自动微分

终于说到自动微分了。

上部分中我们简单聊了聊BP(反向传播，具体可以看CS231n相关的链接)，为什么BP不好呢？是因为其每一步都会保存了上一步中的计算出来的缓冲数据，这样在每次进行反向传播的时候占用的内存比较高。

而拿自动微分来说，自动微分的核心概念是延迟计算的思想。

首先，我们同样选取出一个目标函数，然后求输出值对两个权重参数($W_1、W_2$)的导数。

首先我们求出1/x的导数$-1/x_2$，注意，此时我们仅仅是将这个表达式写了出来，并没有计算这个$-1/x_2$的值。

接下来是前项与1的加法操作的导数，这个导数很显然是1。

再对exp这个op进行求导，同样可以得到导数exp(x)。

最后我们可以得到$J$对$W_1$的偏导，注意这里的偏导是将之前所有得到的式子全部将相应的op的输入代入，然后依次相乘起来。

对于$W_2$来说也同样如此。

可以发现，自动求导和BP最主要的区别是，自动求导没有提前计算出每个算子的导数，仅仅是将计算式子表达出来而已，直到我们确实需要求这个值的时候，整个计算流程才算开始执行。

这是一种延迟计算的思想，我们将所有要计算的路线都规划好之后再进行计算，这样一是可以不用提前计算出中间变量，二是可以根据我们导出的计算节点的拓扑关系进行一些优化之类的工作，总的来说比较灵活。

实现Autodiff算法

那我们尝试实现这样一个autodiff的算法，这也是CSE 599W: Systems for ML这门课中的assignment-1。阅读下文之前可以先看看这个assignment的相关介绍。

首先我们讲一下算法的基本流程：

我们根据目标函数$f = (e^x+1)*e^x$去求每个变量$x_1,x_2,x_3,x_4$的导数。相应的伪代码如下：

这里我们是逆序推导的方式，也就是从输出端一直推导到输入端，而每一个节点都称之为node，例如$x_4=x_2*x_3$，这时$x_4$就相当于一个node，要注意这个node与这个节点的op也有关系，比如这个的op就是mul也就是乘法。

然后我们计算输出的导数(其实就是$x_4$的导数)，显然是1，因为$x_4$就是输出么，这里的node_to_grad储存的是每个node对应的导数。

这里$\overline{x_4}$就表示$x_4$这个node的输出导数。

然后我们走到了第一个节点，$x_4$，我们接着往上走，准备去求$x_4$两个输入的导数($x_2,x_3$)，再提醒一下，这里是逆序。

grad <- sum partial adjoints from output edges这句执行的语句是将目前关于这个node的所有输出的边缘导数相加起来，这样说可能有点难理解。我们之后会说。
然后我们求$x_4$两个输入node的导数，因为是乘法操作，所以$x_2$的导数是$x_3$而$x_3$的导数是$x_2$。

从图中可以看到$x_3$的相对应的导数是$\overline{x_2}^1$(为什么是$\overline{x_2}^1$，后面解释)，$x_2$的相对应的导数是$\overline{x_3}$，而且红色连线的乘号也表示，如果要求$x_2$的导数，则$df/dx_2=\overline{x_4}*\overline{x_3}$。同理，另一个也是如此。

然后我们将已经求出的所有关于node的导数先放到node_to_grad这个map中，以便后续使用。

接下来我们将$x_3$的输出边缘导数都加起来，得到grad，然后用做$x_3$的输入的upstream gradient。

因为$x_3$的输入为$x_2+1$，这里第二次对$x_2$进行求导，与上次不同，这次是在$x_2+1$这个式子中。求出的导数记为$\overline{x_2}^2$。

这里就对$x_2$进行求导后，得到了$x_3=x_2+1$这个式子中，$x_2$的导数。为什么要专门强调一下这个呢？是因为我们的$x_2$在两个node中都出现了，我们拿最开始的目标函数来展示：$f = (e^x+1)*e^x$，其中$x_2$这个node就相当于$e^x$这个op，可以发现$e^x$这个op在这个式子中出现了两次，分别是$(e^x+1)$与$e^x$中，利用求导法则中的乘法法则，分别得到了两个导数值。

这时关于$x^2$的两个导数为$\overline{x_2}^1$和$\overline{x_2}^2$。

将$x^2$的两个导数为$\overline{x_2}^1$和$\overline{x_2}^2$进行求和也就是$\overline{x_2}$，可以看到下面的红色连接的计算路径$\overline{x_2} = \overline{x_4} * \overline{x_2}^1 + \overline{x_2}^1(\overline{x_4}*\overline{x_3})$

对于$x_2$的输入$x_1$，就是本身乘以upstream gradient即可。

将$x_1$的导数结果也添加进去。

最后，由于$x_1$就是输入，整个阶段的自动微分结束。

相关代码

相关的实现代码在assignment-1中，但是需要自己完成重要的部分，这里挑比较重要的简单说一下:

Node是每个计算节点，之后的每个具体的node都要去继承这个类。这个类重载了加法和乘法操作。

lass Node(object):
    """Node in a computation graph."""

    def __init__(self):
        """Constructor, new node is indirectly created by Op object __call__ method.

            Instance variables
            ------------------
            self.inputs: the list of input nodes.
            self.op: the associated op object,
                e.g. add_op object if this node is created by adding two other nodes.
            self.const_attr: the add or multiply constant,
                e.g. self.const_attr=5 if this node is created by x+5.
            self.name: node name for debugging purposes.
        """
        self.inputs = []
        self.op = None
        self.const_attr = None
        self.name = ""

    def __add__(self, other):
        """Adding two nodes return a new node."""
        if isinstance(other, Node):
            new_node = add_op(self, other)
        else:
            # Add by a constant stores the constant in the new node's const_attr field.
            # 'other' argument is a constant
            new_node = add_byconst_op(self, other)
        return new_node

    def __mul__(self, other):
        """Multiplying two nodes return a new node."""
        if isinstance(other, Node):
            new_node = mul_op(self, other)
        else:
            new_node = mul_byconst_op(self, other)
        return new_node

    # Allow left-hand-side add and multiply.
    __radd__ = __add__
    __rmul__ = __mul__

    def __str__(self):
        """Allow print to display node name."""
        return self.name

    __repr__ = __str__

最重要的就是这个自动求导的代码：

def gradients(output_node, node_list):
    """Take gradient of output node with respect to each node in node_list.
    Parameters
    ----------
    output_node: output node that we are taking derivative of.
    node_list: list of nodes that we are taking derivative wrt.
    Returns
    -------
    A list of gradient values, one for each node in node_list respectively.
    """
    # a map from node to a list of gradient contributions from each output node
    # 辅助map 用于存放所有node对应的未整理的grad
    node_to_output_grads_list = {}
    # Special note on initializing gradient of output_node as oneslike_op(output_node):
    # We are really taking a derivative of the scalar reduce_sum(output_node)
    # instead of the vector output_node. But this is the common case for loss function.
    node_to_output_grads_list[output_node] = [oneslike_op(output_node)]
    # a map from node to the gradient of that node
    # 整理后的 对应node的grad map
    node_to_output_grad = {}
    # Traverse graph in reverse topological order given the output_node that we are taking gradient wrt.
    # 我们首先使用dfs先序遍历一遍表达式 之后再逆序一下 就得到反向的元素顺序
    reverse_topo_order = reversed(find_topo_sort([output_node]))

    for node in reverse_topo_order:

        # 一些变量的导数可能有多个部分 需要将计算出来的多项式导数加起来
        grad = sum_node_list(node_to_output_grads_list[node])
        # 最终将计算好的导数存放到 node_to_output_grad 这个 map 中
        node_to_output_grad[node] = grad
        # 根据上一个算子传递过来的 grad 与当前这个 算子node 得到 inputs 的 grads
        input_grads = node.op.gradient(node, grad)

        for id in range(len(node.inputs)):
            if node.inputs[id] not in node_to_output_grads_list:
                node_to_output_grads_list[node.inputs[id]] = [input_grads[id]]
            else:
                node_to_output_grads_list[node.inputs[id]].append(input_grads[id])

    grad_node_list = [node_to_output_grad[node] for node in node_list]
    return grad_node_list

需要注意，只要没有实际执行，那么所有的node都是为计算状态，都只是一系列计算式子 没有实际进行计算。

这里我们简单演示一下，目标函数$y = x_1*x_2+x_1$中，对$x_1$和$x_2$的导数，其中grad_x1和grad_x1为我们需要所求的结果。

x1 = ad.Variable(name="x1")
x2 = ad.Variable(name="x2")

y = x1 * x2 + x1

grad_x1, grad_x2 = ad.gradients(y, [x1, x2])

在debug中可以看到：

因为在深度学习中，最后的损失函数值都是一个数，所以我们求损失对权重的梯度的时候，首先将损失的向量转化为单个的数字，然后用这个数字对权重求梯度，最后得到的值是单个值，不是向量也不是数组。

总结

最后拿两张图来总结下之前说到的几种微分方法，字比较小，建议点开查看：

参考资料

• Clad-DIANA-HEP
• DLSYS-lecture4
• https://blog.csdn.net/qq_40878431/article/details/85942864
• https://medium.com/machine-intelligence-report/how-do-neural-networks-work-57d1ab5337ce
• Automatic Differentiation
  in Machine Learning: a Survey