本系列文章面向程序员，希望通过Image Caption Generation，一个有意思的具体任务，深入浅出地介绍深度学习的知识，涉及到很多深度学习流行的模型，如CNN，RNN/LSTM，Attention等。本文为第三篇。
作者：李理，目前就职于环信，即时通讯云平台和全媒体智能客服平台，在环信从事智能客服和智能机器人相关工作，致力于用深度学习来提高智能机器人的性能。
相关文章：
从Image Caption Generation理解深度学习（part I）
从Image Caption Generation理解深度学习（part II）
2.2.5 反向传播算法的推导
前面我们用很简单的几十行python代码基本上完成了一个多层神经网络。但是还差最重要的部分，那就是计算loss function对参数的偏导数，也就是反向传播算法。下面我们来仔细的完成公式的推导，以及接下来会讲怎么用代码来实现。这一部分数学公式多一些，可能很多读者会希望跳过去，不过我还是建议大家仔细的阅读，其实神经网络用到的数学相比svm,bayes network等机器学习算法，已经非常简单了。请读者阅读的时候最好准备一支笔和几张白纸，每一个公式都能推导一下。如果坚持下来，你会觉得其实挺简单的。
(1) feedforward阶段的矩阵参数表示和计算
之前我们讨论的是一个神经元的计算，而在代码里用到的却是矩阵向量乘法。而且细心的读者会发现我们在构造参数矩阵weights的时候，行数和列数分别是后一层的节点数和前一层的节点数。这似乎有点不自然，为什么不反过来呢？看过下面这一部分就会明白了。
首先我们熟悉一下第L（因为小写的L和1太像，所以我用大写的L）层的参数w_jk。它表示第L-1层的第k个神经元到第L层的第j个神经元的权重。比如第3层的w_24，参考上面的图，它表示的是第2层的第4个神经元到第3层的第二个神经元。
对bias和激活函数后的结果a也采用类似的记号，如下图所示。
b_32表示第2层的第3个神经元的bias，而a_13第3层的第1个神经元的激活。
使用上面的记号，我们就可以计算第L层的第j个神经元的输出a_jl：
第L层的第j个神经元的输入是L-1层的a_1,a_2,...；对应的权值是w_j1,w_j2,...；bias是b_jL。所以a_jL就是上面的公式，k的范围是从1到第L-1层的神经元的个数。
为了用矩阵向量乘法来一次计算第L层的所有神经元的输出，我们需要定义第L层的参数矩阵w_l，它的大小是m*n，其中m是第L层的神经元个数；而n则是第L-1层的个数。它的第i行第j列就是我们上面定义的w_jk。此外我们还要定义向量b_l，它的大小是m（也就是第L层神经元的个数），它的第j个元素就是我们上面定义的b_j。
最后，我们定义element-wise的函数，比如f(x) = x^2，如果输入是一个向量，那么结果是和输入一样大小的向量，它的每个元素是对输入向量的每一个元素应用这个函数的结果。
有了上面的定义，我们就可以一次计算出第L层的输出（一个长度为m的向量）
下面是对上面这个公式的详细证明（说明）：
我们需要证明的是向量aL的第j个元素就是前面的a_jL
此外，为了方便后面的求解，我们把加权累加和也用一个符号z_l来表示。
其中，它的第j个元素就是第L层的第j个神经元的加权累加和：
这样a_l就可以简单的对z_l的每个元素计算激活函数
现在我们再回顾一下feedforward的代码就非常直观了：
def feedforward(self, a):"""Return the output of the network if a is input."""for b, w in zip(self.biases, self.weights): a = sigmoid(np.dot(w, a)+b)return a
传给函数feedforward的参数a就是输入向量x，第一层就是x，第二层就是第一个隐层，每一层的计算就是非常简单的参数矩阵w_l乘以上一层的激活a_l-1在加上b_l，然后用激活函数计算。
初始化的时候w的大小是 (后一层的神经元个数) * (前一层的神经元个数），再回顾一下初始化参数的代码：
# sizes = [784, 30, 10]def __init__(self, sizes):self.num_layers = len(sizes) self.sizes = sizes self.biases = [np.random.randn(y, 1) for y in sizes[1:]] self.weights = [np.random.randn(y, x)for x, y in zip(sizes[:-1], sizes[1:])]
x, y in zip(sizes[:-1], sizes[1:]) x是第一层到最后倒数第二层，y是第二层到最后一层，比如上面的sizes=[784, 30, 10]
x是[784, 30], y是[30, 10]，注意随机的矩阵是(y,x)，所以self.weights是两个矩阵，大小分别是30*784和10*30
(2) 关于损失函数C的两个假设
1. 损失函数是每个训练数据的损失的平均
也就是C是这样的形式：
对于之前我们使用的MSE损失函数，这是满足的。我们使用batch的梯度下降的时候需要求C对参数w的偏导数，因为损失函数是每个训练数据的损失的平均，所以我们只需要求每个数据的偏导数，然后加起来平均就行。这个假设几乎所有的损失函数都是满足的【我是没见过损失函数不满足这个条件】
损失函数是最后一层输出的函数
这个条件几乎常见的损失函数都是这样的，我们之前时候的MSE就是计算最后一层的输出aL和正确的y(one-hot)的均方误差，显然是满足的。
(3) Hadamard product
这个名字看起来很复杂，其实很简单，就是两个向量elementwise的乘法。看一个例子就清楚了：
(4) 反向传播算法(back propagation)的4个公式
回顾一下，我们之前说了，梯度下降其实最核心的问题就是求损失函数对每一个参数的偏导数。那我们就直接一个一个求好了，为什么又要搞出一个反向传播算法呢？其实这个算法在不同的领域被不同的人重复“发现”过很多次，有过很多不同的名字，最本质的应该就是逆向求导（reverse-mode differentiation）或者叫做自动求导（automatic differentiation）。自动求导（AD）是非常通用的一种求偏导数的方法，很早就在流体力学和大气物理等领域使用，反向传播算法可以认为是AD在神经网络中的应用。不过最早发现这个算法的人（是谁最早好像还有点争议）并不是先知道AD可以直接用于神经网络，他发现这个算法是基于错误的反向传播而得到的，所有命名为（错误的）反向传播算法。后面我们会讲到AD，这是一个强大的算法，任何一个函数，你能把它分解成有向无环图的计算图【函数一般都能分解成一些无依赖的最基础的变量的复合函数，因此肯定可以表示成这样一个有向无环图】，然后每个节点都表示一个函数。只要你能求出这个函数在特定点的梯度【也就是这个函数对所以自变量的偏导数】（不需要求解析的偏导数，当然很多情况，这些函数都是能直接求出解析解，然后代入这个特定点就行，但理论上我们是可以用其他方法，比如数值梯度近似来求的），就能自动的计算损失函数对每一个参数的偏导数（也是在这个点的），而且只要反向根据拓扑排序遍历这个图一次就行，非常高效和简单。后面我们会详细的介绍AD。这个方法非常通用，TensorFlow的核心就是AD。使用AD的框架就比较灵活，我想“创造”一种新的网络结构，我又不想【其实更可能是不会】推导出梯度的公式，那么我只需要把我的网络能用这样一个有向无环图表示就行。当然节点必须要能够求出梯度来，一般我们的函数比如矩阵的运算，卷积等等TensorFlow都封装好了——它把它叫做一个op。我们只需要搭积木一样把这个计算图定义出来，TensorFlow就自动的能根据AD计算出损失函数对所有参数的梯度来了。当然如果你要用到一个TensorFlow没有的op，那你就需要根据它的规范实现这个op，一个op最核心的接口就是两个，一个是输入x，求f(x)；另一个就是求f在某个x0点的梯度。
不过这里，我们还是沿着神经网络的发展历史，从错误的反向传播角度来理解和推导这个算法。
首先，我们会对每一个神经元比如第L层的第j个，都定义一个错误δ_jL
也就是损失函数对z也就是线性累加和的偏导数。为什么定义这样一个东西呢？我们假设在第L层的第j个神经元上有一个精灵（Daemon）
当这个神经元得到来自上一次的输入累加计算出z_jL的时候，它会恶作剧的给一点很小的干扰Δz_jL。原来它应该输出的是σ(z_jL)，现在变成了σ(z_jL +Δz_jL)。这个微小的变化逐层传播，最终导致损失函数C也发生如下的变化：
这个其实就是导数的直觉定义：微小的Δx引起微小的Δy，Δy/Δx约等于导数。
不过这个精灵是个好精灵，它想帮助我们减少损失。 当
大于0的时候，它让Δz_jL小于0，反之当它小于0的时候它让Δz_jL大于0。这样
总是小于0
因此我们的loss就会变小。而其绝对值越大，我们的损失减少的越多。
当然你会说为什么不能让Δz_jL非常大，这样我们的损失总是减少很多？可惜这个精灵是个数学家，它说如果Δx太大，那么Δy=df/dx *Δx就不准确了。
所以我们可以这样认为：它就是第L层的第j个神经元“引起”的“错误”。如果绝对值大，则它的“责任”也大，它就得多做出一些调整；反之如果它趋近于0，说明它没有什么“责任”，也就不需要做出什么改变。
因此通过上面的启发，我们定义出δ_jL来。
接下来我们逐个介绍反向传播算法的4个公式。
公式1. 第L层（最后一层） 的错误
这个公式的第一项，就是损失C对a_jL的导数，它越大，说明C受a_jL的影响也就越大，如果有了错误，第a_jL的“责任”也就越大，错误也就越大。第二项是a_jL受z_jL的影响。两者乘起来就是z_jL对最终损失的影响，也就是它的“责任”的大小。
这个公式很好计算，首先第二项就是把z_jL的值（这个在feedforward节点就算出来并存储下来了）代入σ'(x)。如果σ是sigmoid函数，我们前面也推导过它的导数：σ’(x)=σ(x)*(1-σ(x))。第一项当然依赖于损失函数的定义，一般也很好求。比如我们的MSE损失：
具体的推导我在纸上写了一下，虽然很简单，我们也可以练练手，尤其是对于求和公式的展开，希望大家能熟悉它，以后的推导我可能就不展开求和公式了，你需要知道求和公式里哪些项是和外面的自变量无关的。
公式BP1是elementwise的，我们需要变量j来计算每一个δ_jL。我们也可以把它写成向量的形式，以方便利用线性代数库，它们可以一次计算向量或者矩阵，可以用很多技术利用硬件特性来优化（包括GPU，SSE等）速度。
右边δ'(z_L)很容易理解，左边的记号可能有些费解，其实我们把∇aC当成一个整体就好了，它是一个向量，第一个元素是∂C/∂a_1L，第二个就是∂C/∂a_2L，…
如果算上函数C是MSE的话，上面的公式就可以简化成：
公式2. 第l层（非最后一层） 的错误
等下我们会证明这个公式，不过首先我们来熟悉一下公式。如果我们想“背”下这个公式的话，似乎看起来比第一个BP1要复杂很多 。我们先检查一下矩阵和向量的维度，假设l+1层有m个元素，l层n个。则w_l+1的大小是m*n，转置之后是n*m, δ_l+1的大小是n*1，所以矩阵相乘后是m*1，这和δ_l是一样的，没有问题。
接下来我们仔细观察一下BP2这个公式，首先第二项σ'(z_l)和前面的含义一样，代表a_l对于z_l的变化率。
而第一项复杂一点，我们知道第l层的第j个神经元会影响第l+1层的所有神经元，从而也影响最终的损失C。这个公式直接给了一个矩阵向量的形式，看起来不清楚，所以我在草稿纸上展开了：
最终第L层的第j个神经元的损失就是如下公式：
这下应该就比较清楚了，第l层的第j个神经元的损失，就是把l+1层的损失“反向传播”回来，当然要带上权重，权重越大，“责任”也就越大。
如果要“背”出这个公式也没有那么复杂了，先不看σ'(z_l)，第一项应该是矩阵w_l+1乘以δ_l+1。由于矩阵是m*n，而
向量δ_l+1是m*1，为了能让矩阵乘法成立，那么就只能把w转置一下，变成n*m，然后就很容易记住这个公式了。
注意，BP2的计算是从后往前的，首先根据BP1，最后一层的δ_L我们已经算出来了，因此可以向前计算L-1层的δ_L-1，
有了δ_L-1就能计算δ_L-2，…，最终能算出第一个隐层（也就是第2层）δ_1来。
公式3. 损失函数对偏置b的梯度
这前面费了大力气求δ_l，不要忘了我们的最终目标是求损失函数对参数w和b的偏导数，而不是求对中间变量z的偏导数。
因此这个公式就是对b的偏导数。
或者写成向量的形式：
∂C/∂b就是δ！
公式4. 损失函数对w的梯度
或者参考下图写成好记的形式：
也就是说对于一条边w_jkL，∂C/∂w_ij就是这条边射出的点的错误δ乘以进入点的激活。非常好记。
我们把这四个公式再总结一下：
(5) 这四个公式的证明
首先是BP1，请参考下图：
然后是BP2：
这里用到了chain rule，其实也非常简单和直观，就是复合函数层层组合。最简单的方法就是用图画出来，比如y最终
是x的函数，我们要求∂y/∂x，如果y是u,v的函数，然后u,v才是x的函数，那么我们把变量x,y,u,v都画成图上的点，y是u，v的函数，那么我们画上从u和v到y的边，同样，我们画上从x到u和v的边，然后从y到x的每条路径，我们经过的边都是一个偏导数，我们把它累加起来就行【这其实就是后面我们会讲的AD】。因此∂y/∂x=∂y/∂u * ∂u/∂x +∂y/∂v * ∂v/∂x。
剩下的BP3和BP4也非常类似，我就不证明了。
反向传播算法
1. a_1 = 输入向量x
2. Feedforward 根据公式
和
计算z_l和a_l并存储下来（反向传播时要用的）
3. 计算最后一层的错误
向前计算所有层的错误
计算损失对所有参数的偏导数
2.2.6 代码实现反向传播算法
我们已经把公式推导出来了，那怎么用代码实现呢？我们先把代码复制一下，然后说明部分都是作为代码的注释了，
请仔细阅读。
class Network(object): def update_mini_batch(self, mini_batch, eta): # mini_batch是batch大小，eta是learning rate nabla_b = [np.zeros(b.shape) for b in self.biases] # 构造和self.biases一样大小的向量，比如前面的例子 sizes=[784,30,10]，则 # nabla_b是两个向量，大小分别是30和10 nabla_w = [np.zeros(w.shape) for w in self.weights] # 构造和self.weights一样大小的矩阵，比如前面的例子 sizes=[784,30,10]，则 # nabla_w是两个矩阵，大小分别是30*784和10*30 for x, y in mini_batch: #对于每个训练样本x和y delta_nabla_b, delta_nabla_w = self.backprop(x, y) # 用backprop函数计算损失函数对每一个参数的偏导数。 # backprop函数下面会详细讲解 nabla_b = [nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)] # 把返回的对b偏导数累加到nabla_b中 nabla_w = [nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)] # 把返回的对w的偏导数累加到nabla_w中 self.weights = [w-(eta/len(mini_batch))*nw for w, nw in zip(self.weights, nabla_w)] # 计算完一个batch后更新参数w self.biases = [b-(eta/len(mini_batch))*nb for b, nb in zip(self.biases, nabla_b)] # 更新b... def backprop(self, x, y): # 输入是x和y，返回损失函数C对每个参数w和b的偏导数 # 返回的格式是两个元组，第一个是b的偏导数，第二个是w的。 nabla_b = [np.zeros(b.shape) for b in self.biases] # 构造和self.biases一样大小的向量，比如前面的例子 sizes=[784,30,10]，则 # nabla_b是两个向量，大小分别是30和10 nabla_w = [np.zeros(w.shape) for w in self.weights] # 构造和self.weights一样大小的矩阵，比如前面的例子 sizes=[784,30,10]，则 # nabla_w是两个矩阵，大小分别是30*784和10*30 # feedforward activation = x activations = [x] # 用一个list保存所有层的激活，下面backward会有用的 zs = [] # 同样的用一个list保存所有层的加权累加和z，下面也会用到。 #下面这段代码在feedward也有，不过那里是用来predict用的不需要保存zs和activations for b, w in zip(self.biases, self.weights): z = np.dot(w, activation)+b zs.append(z) activation = sigmoid(z) activations.append(activation) # backward pass #1. 首先计算最后一层的错误delta，根据公式BP1，它是损失函数对a_L的梯度乘以σ'(z_L) # sigmoid_prime就是σ'(z_L)，而∂C/∂a_L就是函数cost_derivative，对于MSE的损失函数， # 它就是最后一层的激活activations[-1] - y delta = self.cost_derivative(activations[-1], y) * \ sigmoid_prime(zs[-1]) # 2. 根据公式BP3,损失对b的偏导数就是delta nabla_b[-1] = delta # 3. 根据公式BP4,损失对w的偏导数时delta_out * activation_in # 注意，我们的公式BP4是elementwise的，我们需要写成矩阵向量的形式 # 那怎么写呢？我们只需要关心矩阵的大小就行了。 # 假设最后一层有m(10)个神经元，前一层有n(30)个， # 则delta是10*1, 倒数第二层的激活activations[-2]是30*1 # 我们想求的最后一层的参数nabla_w[-1]是10*30，那么为了能够正确的矩阵乘法， # 只要一种可能就是 delta 乘以 activations[-2]的转置，其实也就是向量delta和activations[-2]的外积 nabla_w[-1] = np.dot(delta, activations[-2].transpose()) # 接下来从倒数第二层一直往前计算delta，同时也把对w和b的偏导数求出来。 # 这里用到一个比较小的trick就是python的下标是支持负数的，-1表示最后一个元素，-2是倒数第二个 # l表示倒数第l层，2就表示倒数第2层，num_layers - 1就表示顺数第2层（也就是第1个隐层） # 比如我们的例子：sizes=[784, 30, 10]，那么l就是从2到3（不包含3），l就只能是2，页就是第1个（也是唯一的一 # 个）隐层 for l in xrange(2, self.num_layers): # 倒数第l层的z z = zs[-l] # 计算σ'(z_l) sp = sigmoid_prime(z) # 根据BP2，计算delta_l,注意weights[-l+1]表示倒数第l层的下一层 delta = np.dot(self.weights[-l+1].transpose(), delta) * sp # 同上，根据BP3 nabla_b[-l] = delta # BP4，矩阵乘法参考前面的说明 nabla_w[-l] = np.dot(delta, activations[-l-1].transpose()) return (nabla_b, nabla_w)
2.2.7 为什么反向传播算法是一个高效的算法？
分析完代码，我们发现一次backprop函数调用需要feedforward一次，网络有多少边，就有多少次乘法，有多少个点就有多少次加分和激活函数计算（不算第一层输入层）。反向计算也是一样，不过是从后往前。也就是说这是时间复杂度为O(n)的算法。
如果我们不用反向传播算法，假设我们用梯度的定义计算数值梯度。对于每一个参数wj，
我们都用公式 limit (f(w1, w2, …, wj+Δ wj, …) - f(w1, w2, …, wj, …)/Δwj
f(w1, w2, wj, …)只需要feedforward一次，但是对于每个参数wj，都需要feedforward一层来计算f(w1, w2, …, wj+Δ wj, …)，它的时间复杂度是O(n)，那么对所有的参数的计算需要O(n^2)的时间复杂度。
假设神经网络有1百万个参数，那么每次需要10^12这个数量级的运算，而反向传播算法只需要10^6，因此这个方法比反向传播算法要慢1百万倍。
