这里把按 [1] 推导的BP算法（Backpropagation）步骤整理一下。突然想整理这个的原因是知乎上看到了一个帅呆了的求矩阵微分的方法（也就是 [2]），不得不感叹作者的功力。[1] 中直接使用矩阵微分的记号进行推导，整个过程十分简洁。而且这种矩阵形式有一个非常大的优势就是对照其进行编程实现时非常方便。
但其实用标量计算推导也有一定的好处，比如可以清楚地知道某个权重是被谁所影响的。
前向传播过程：多层Logistic回归
记号约定：
$L$：神经网络的层数。输入层不算。
$n^l$：第 $l$ 层神经元的个数。偏置神经元不算在内。
$W^{l}\in\mathbb R^{n^l\times n^{l-1}}$：第 $l-1$ 层到第 $l$ 层的权重矩阵。其中，$W_{ij}^{(l)}$ 表示第 $l-1$ 层第 $j$ 个神经元到第 $l$ 层第 $i$ 个神经元的权重。
$\textbf b^{(l)}\in\mathbb R^{n^l}$ ：第 $l$ 层的偏置。
$\textbf z^{(l)}\in\mathbb R^{n^l}$ ：第 $l$ 层各个神经元的输入。
$f_l(\cdot)$ ：第 $l$ 层的激活函数。对于分类任务来说，最后一层为softmax。
$\textbf a^{(l)}\in\mathbb R^{n^l}$ ：第 $l$ 层各个神经元的输出，也就是活性值。对于输入层（第0层），$\textbf a^{(0)}=\textbf x$。
下图给出了相应的图示。
图片来源：[1]
下面可以开始推导了。首先是前向传播，计算网络输出：
$$\textbf z^{(l)}=W^{(l)}\textbf a^{(l-1)}+\textbf b^{(l)}$$
$$\textbf a^{(l)}=f_l(\textbf z^{(l)})=f_l(W^{(l)}\textbf a^{(l-1)}+\textbf b^{(l)})$$
$$\textbf a^{(L)}=f(\textbf x;W,\textbf b)$$
整个网络实际上就可以看作是一个函数 $f$ ，1989年已经有人证明了单隐层网络可以逼近任意函数。
既然是完成从输入到输出的映射，那么网络除了可以看作分类器之外，还可以看作是特征提取器，再将网络的输出作为分类器的输入，相当于把原特征 $\textbf x$ 映射成了新的特征 $\textbf a^{(L)}$ 。
如果把网络看作分类器，那么最后一层的神经元个数 $n^L$ 应等于类别个数 $C$ ，且 $f_L(\cdot)=\text{softmax}(\cdot)$ ，归一化成概率分布：
$$\hat{\textbf y}=\textbf a^{(L)}=\text{softmax}(\textbf z^{(L)})=\frac{\exp(\textbf z^{(L)})}{\textbf 1^{\top}\exp(\textbf z^{(L)})}\in\mathbb R^C$$
需要注意的是分子是列向量，分母是标量。$\hat{\textbf y}$ 的每一维 $\hat{\textbf y}_i$ 表示的是网络 $f$ 给出的样本 $\textbf x$ 属于第 $i$ 类的概率。
给定样本 $(\textbf x,\textbf y)$ ，其中 $\textbf y\in\mathbb R^C$ 是one-hot向量，那么使用交叉熵损失函数，得到网络对一个样本的损失为
$$\mathcal L(\textbf y,f(\textbf x;W,\textbf b))=-\textbf y^{\top}\ln\hat{\textbf y}$$
设训练样本数为 $N$ ，那么经验风险为
$$R=\frac1N\sum_{i=1}^N\mathcal L(\textbf y_i,f(\textbf x_i;W,\textbf b))$$
如果使用权重矩阵的Frobenius范数作为正则项
$$||W||_F=\biggl(\sum_{l=1}^L\sum_{i=1}^{n^l}\sum_{j=1}^{n^{l-1}}(W_{ij}^{(l)})^2\biggr)^{\frac12}$$
那么结构风险为
$$R=\frac1N\sum_{i=1}^N\mathcal L(\textbf y_i,f(\textbf x_i;W,\textbf b))+\frac12\lambda||W||_F^2$$
一般来说，偏置是不加正则的。
现在考虑结构风险最小化，使用梯度下降法来更新权重矩阵和偏置。只需要求取模型对一个样本的损失的梯度 $\dfrac{\partial \mathcal L(\textbf y,f(\textbf x;W,\textbf b))}{\partial W^{(l)}}$（因为对多个样本的梯度就是对一个样本的梯度的累加）：
$$W^{(l)}=W^{(l)}-\alpha \frac{\partial R}{\partial W^{(l)}}$$
$$\frac{\partial R}{\partial W^{(l)}}=\frac1N\sum_{i=1}^N\frac{\partial \mathcal L(\textbf y,f(\textbf x;W,\textbf b))}{\partial W^{(l)}}+\lambda W^{(l)}$$
对偏置的就不写了。
反向传播（Backpropagation，BP）算法
好了，墨迹了半天，重点终于来了。下面就是要求解 $\dfrac{\partial\mathcal L(\textbf y,f(\textbf x;W,\textbf b))}{\partial W^{(l)}}$ 、$\dfrac{\partial\mathcal L(\textbf y,f(\textbf x;W,\textbf b))}{\partial \textbf b^{(l)}}$ 。为了简单起见，直接记作 $\dfrac{\partial \mathcal L}{\partial W^{(l)}}$ 、$\dfrac{\partial \mathcal L}{\partial \textbf b^{(l)}}$ 。
我所有的符号都是沿用的 [1] ，后面的推导过程也是。
首先摆上一些公式复习一下：
$$\frac{\partial A^{\top}\textbf x}{\partial\textbf x}=\frac{\partial \textbf x^{\top}A}{\partial\textbf x}=A$$
$$\frac{\partial \textbf y^{\top}\textbf z}{\partial\textbf x}=\frac{\partial \textbf y}{\partial\textbf x}\textbf z+\frac{\partial \textbf z}{\partial\textbf x}\textbf y$$
$$\frac{\partial \textbf y^{\top}A\textbf z}{\partial\textbf x}=\frac{\partial \textbf y}{\partial\textbf x}A\textbf z+\frac{\partial \textbf z}{\partial\textbf x}A^{\top}\textbf y$$
$$\frac{\partial y\textbf z}{\partial\textbf x}=\frac{\partial y}{\partial\textbf x}\textbf z^{\top}+y\frac{\partial \textbf z}{\partial\textbf x}$$
$$\frac{\partial \text{tr}AB}{\partial A}=B^{\top}\quad\quad\frac{\partial \text{tr}AB}{\partial A^{\top}}=B$$
$$\frac{\partial f(A)}{\partial A^{\top}}=(\frac{\partial f(A)}{\partial A})^{\top}$$
然后是链式法则（chain rule）：
$$\frac{\partial \textbf z}{\partial \textbf x}=\frac{\partial \textbf y}{\partial \textbf x}\frac{\partial \textbf z}{\partial \textbf y}$$
$$\frac{\partial z}{\partial X_{ij}}=(\frac{\partial z}{\partial\textbf y})^{\top}\frac{\partial\textbf y}{\partial X_{ij}}$$
$$\frac{\partial z}{\partial X_{ij}}=\text{tr}\biggl((\frac{\partial z}{\partial Y})^{\top}\frac{\partial Y}{\partial X_{ij}}\biggr)$$
以及逐元素计算的函数 $f$ （其导函数为 $f'$ ）及其梯度：
$$\frac{\partial f(\textbf x)}{\partial\textbf x}=\text{diag}(f'(\textbf x))$$
这里需要声明一点，上面这套计算体系里，都是使用的分母布局，也就是说：维度为 $p$ 的列向量对维度为 $q$ 的列向量求导后得到的矩阵维数为 $q\times p$ 。关于矩阵求导，可以参考 [3] 。
又墨迹了半天，下面开始求 $\dfrac{\partial \mathcal L}{\partial W^{(l)}}$ 、$\dfrac{\partial \mathcal L}{\partial \textbf b^{(l)}}$ 。
BP的四个重要公式
下面这第一步在我认为是最重要的一步，因为它的思路很巧妙 —— 现在的目标是求标量对矩阵的导数，直接用chain rule拆成两块的话会发现：一块是标量对向量的导数，好办；还有一块是向量对矩阵的导数，不好办。所以 [1] 使用了一种很巧的方式，求标量对矩阵单个元素（标量）的导数，然后再归纳成矩阵形式。
根据链式法则，
$$\frac{\partial\mathcal L}{\partial W_{ij}^{(l)}}=(\frac{\partial\mathcal L}{\partial \textbf z^{(l)}})^{\top}\frac{\partial \textbf z^{(l)}}{\partial W_{ij}^{(l)}}$$
现在定义误差项 $\delta^{(l)}=\dfrac{\partial\mathcal L}{\partial \textbf z^{(l)}}\in\mathbb R^{n^l}$ ，用来表征第 $l$ 层的神经元对于误差的敏感程度。
先求 $\dfrac{\partial \textbf z^{(l)}}{\partial W_{ij}^{(l)}}$ ：
$$\frac{\partial \textbf z^{(l)}}{\partial W_{ij}^{(l)}}=\frac{\partial W^{(l)}\textbf a^{(l-1)}}{\partial W_{ij}^{(l)}}=(0,...,a_j^{(l-1)},...,0)^{\top}$$
其中只有第 $i$ 行的元素非零。既然这样，那么
$$\frac{\partial\mathcal L}{\partial W_{ij}^{(l)}}=\delta_i^{(l)}a_j^{(l-1)}$$
从这个式子可以看出，如果神经元的输出值 $a_j^{(l-1)}$ 很小，那么它连向下一层的权重 $W_{ij}^{(l)}$ 的更新就会很缓慢。
写成矩阵的形式，就是下面的式子：
$$\frac{\partial\mathcal L}{\partial W^{(l)}}=\delta^{(l)}(\textbf a^{(l-1)})^{\top}$$
$$\frac{\partial\mathcal L}{\partial \textbf b^{(l)}}=\delta^{(l)}$$
下面的问题就是如何计算误差项 $\delta^{(l)}$ 。对于误差项的计算，需要分两种情况考虑，一种情况是输出层，另一种情况是隐层。
对于输出层来说，
$$\begin{aligned}\delta^{(L)}&=\frac{\partial \mathcal L}{\partial \textbf z^{(L)}}\\&=\frac{\partial \textbf a^{(L)}}{\partial \textbf z^{(L)}}\frac{\partial \mathcal L}{\partial \textbf a^{(L)}}\\&=\text{diag}(f_L'(\textbf z^{(L)}))\frac{\partial \mathcal L}{\partial \textbf a^{(L)}}\\&=f_L'(\textbf z^{(L)})\odot \frac{\partial \mathcal L}{\partial \textbf a^{(L)}}\end{aligned}$$
对于隐层来说，
$$\begin{aligned}\delta^{(l)}&=\frac{\partial \mathcal L}{\partial \textbf z^{(l)}}\\&=\frac{\partial \textbf a^{(l)}}{\partial \textbf z^{(l)}}\frac{\partial\textbf z^{(l+1)}}{\partial \textbf a^{(l)}}\frac{\partial \mathcal L}{\partial \textbf z^{(l+1)}}\\&=\text{diag}(f_l'(\textbf z^{(l)}))(W^{(l+1)})^{\top}\delta^{(l+1)}\\&=f_l'(\textbf z^{(l)})\odot \bigl((W^{(l+1)})^{\top}\delta^{(l+1)}\bigr)\end{aligned}$$
从这个式子就可以看出，误差项可以逐层往输入的方向传递，这就是所谓的误差反向传播。
多提一句，这里出现了激活函数的导函数。有以下两个常用关系：
$$\sigma'(\cdot)=\sigma(\cdot)\odot (\textbf 1-\sigma(\cdot))$$
$$\text{softmax}'(\cdot)=\text{softmax}(\cdot)\odot (\textbf 1-\text{softmax}(\cdot))$$
其中 $\sigma(\cdot)$ 代表logistic函数。
所以，BP算法归结到最后其实只有四个式子比较重要：
$$\delta^{(L)}=f_L'(\textbf z^{(L)})\odot \frac{\partial \mathcal L}{\partial \textbf a^{(L)}}$$
$$\delta^{(l)}=f_l'(\textbf z^{(l)})\odot \bigl((W^{(l+1)})^{\top}\delta^{(l+1)}\bigr)$$
$$\frac{\partial \mathcal L}{\partial W^{(l)}}=\delta^{(l)}(\textbf a^{(l-1)})^{\top}$$
$$\frac{\partial \mathcal L}{\partial \textbf b^{(l)}}=\delta^{(l)}$$
BP算法与Delta规则的比较
1. 与二类感知器（只有一层，所以不需要看第二个式子）的Delta规则的比较：BP算法 $\delta^{(L)}$ 表达式中的 $f_L'(\textbf z^{(L)})$ 这项是额外多出来的。这也不难理解，因为二类感知器的最后是经过一个与0.5比较大小的阈值函数，这里则是Softmax函数（多分类）或Logistic函数（二分类）。
二类感知器中准则函数对参数的梯度就是用误差（被错误分类的样本为-1，被正确分类的样本为0，对应到这里就是 $\dfrac{\partial \mathcal L}{\partial \textbf a^{(L)}}$ ）乘以样本的特征（对应到这里就是上一层输出 $\textbf a^{(L-1)}$ ）。
2. 另外可以看出，网络的各个层的误差项之间是相乘的关系，那么当层数很深的时候就可能会出现梯度消失／梯度爆炸的问题。
我也用标量的形式推导过，过程如下……
无论标量形式还是矩阵形式，关键就是一点 —— chain rule。
输出层的误差项 $\delta^{(L)}$求解
上面四个式子只是形式化地给出了反向传播算法的梯度形式。下面求取：对于分类问题，也就是 $f_L(\cdot)=\text{softmax}(\cdot)$ ，并且使用交叉熵损失函数 $\mathcal L(\textbf y,f(\textbf x;W,\textbf b))=-\textbf y^{\top}\log\hat{\textbf y}$ 时，输出层的误差项 $\delta^{(L)}=\dfrac{\partial \mathcal L}{\partial \textbf z^{(L)}}$ 到底是个什么形式。
其实到这里问题已经很清楚了：这和Softmax回归模型是一样一样的，这里的 $\delta^{(L)}$ 在Softmax回归模型的求梯度过程中是一个中间结果。有三种方法：
（1）普通方法，按部就班地一步步推，可以参考我很早之前写的一篇关于word2vec的博客；
（2）[1] 中的方法（第三章，在Softmax回归那部分）；
（3）[2] 中的方法。
下面是基于 [2] 所给方法的完整求解过程。为了简洁起见，所有的 $\textbf z^{(L)}$ 都记作 $\textbf z$：
[2]介绍的是这样形式的求导：已知矩阵 $X$ ，函数 $f(X)$ 的函数值为标量，求 $\dfrac{\partial f}{\partial X}$ 。一种典型的例子就是求取损失对权重矩阵的导数。
对于一元微积分，$\text{d}f=f'(x)\text{d}x$ ；多元微积分，$\text{d}f=\sum_i\dfrac{\partial f}{\partial x_i}\text{d}x_i=(\dfrac{\partial f}{\partial \textbf x})^{\top}\text{d}\textbf x$；由此建立矩阵导数和微分的联系：
$$\text{d}f=\sum_{i,j}\frac{\partial f}{\partial X_{ij}}\text{d}X_{ij}=\text{tr}((\frac{\partial f}{\partial X})^{\top}\text{d}X)$$
上式第二个等号成立是因为对于两个同阶方阵有 $\text{tr}(A^{\top}B)=\sum_{i,j}A_{ij}B_{ij}$ 。求解的流程就是，先求微分 $\text{d}f$ 表达式，然后再套上迹（因为标量的迹等于标量本身），然后再把表达式 $\text{tr}(\text{d}f)$ 和 $\text{tr}((\dfrac{\partial f}{\partial X})^{\top}\text{d}X)$ 进行比对，进而把 $\dfrac{\partial f}{\partial X}$ 给“挖”出来。
所以，问题就从求梯度转化成了求微分。求微分当然少不了很多法则和技巧，下面随着讲随着介绍。
首先求取 $\text{d} \mathcal L$ 。
$$\begin{aligned} \mathcal L&=-\textbf y^{\top}\ln\frac{\exp(\textbf z)}{\textbf 1^{\top}\exp(\textbf z)}\\&=-\textbf y^{\top}(\textbf z-\ln\begin{pmatrix}\textbf 1^{\top}\exp(\textbf z) \\ \textbf 1^{\top}\exp(\textbf z) \\ \vdots \\ \textbf 1^{\top}\exp(\textbf z)\end{pmatrix})\quad \textbf 1^{\top}\exp(\textbf z)\text{是标量}\\&=\ln(\textbf 1^{\top}\exp(\textbf z))-\textbf y^{\top}\textbf z\end{aligned}$$
根据法则 $\text{d}(g(X))=g'(X)\odot\text{d}X$ 、$\text{d}(XY)=(\text{d}X)Y+X(\text{d}Y)$，可得
$$\text{d}(\ln(\textbf 1^{\top}\exp(\textbf z)))=\frac{1}{\textbf 1^{\top}\exp(\textbf z)}\odot\text{d}(\textbf 1^{\top}\exp(\textbf z))$$
$$\text{d}(\textbf 1^{\top}\exp(\textbf z))=\textbf 1^{\top}\text{d}(\exp(\textbf z))=\textbf 1^{\top}(\exp(\textbf z)\odot\text{d}\textbf z)$$
所以
$$\text{d} \mathcal L=\frac{\textbf 1^{\top}(\exp(\textbf z)\odot\text{d}\textbf z)}{\textbf 1^{\top}\exp(\textbf z)}-\textbf y^{\top}\text{d}\textbf z$$
现在可以套上迹，根据恒等式 $\text{tr}(A^{\top}(B\odot C))=\text{tr}((A\odot B)^{\top}C)=\sum_{i,j}A_{ij}B_{ij}C_{ij}$ ，可得
$$\begin{aligned}\text{d} \mathcal L&=\text{tr}(\frac{(\textbf 1\odot \exp(\textbf z))^{\top}\text{d}\textbf z}{\textbf 1^{\top}\exp(\textbf z)})-\text{tr}(\textbf y^{\top}\text{d}\textbf z)\\&=\text{tr}(\biggl(\frac{(\exp(\textbf z))^{\top}}{\textbf 1^{\top}\exp(\textbf z)}-\textbf y^{\top}\biggr)\text{d}\textbf z)\\&=\text{tr}((\hat{\textbf y}-\textbf y)^{\top}\text{d}\textbf z)\\&=\text{tr}((\frac{\partial\mathcal L}{\partial\textbf z})^{\top}\text{d}\textbf z)\end{aligned}$$
这样就得到了 $\delta^{(L)}$ 的形式：
$$\delta^{(L)}=\frac{\partial \mathcal L}{\partial \textbf z^{(L)}}=\hat{\textbf y}-\textbf y$$
这样就不难看出为什么将 $\delta^{(l)}$ 称为误差项了。
学习率自适应
众所周知，现在有很多算法用来改进SGD，以求更好的训练效果（时间上 & 精度上）。下面这张图总结了包括加动量、自适应学习率等众多改进算法的更新公式，如果想详细学习其中的理论，可以看 [1] ，也可以看Yoshua Bengio的那本最新出版的《Deep Learning》。对于这些算法，我认为没有最好，只有最合适。在自己的任务上何者可以取得最满意的效果，还是要多跑几组实验来确定。反正在TensorFlow里就是改个函数的事情，但是挑出最适合自己任务的那一个，还是挺不容易的。另外，诸如什么“先用 Adam 再用 SGD ”之类的我是没有试过（我都是直接用Adam，逃。。。
图片来源：http://weibo.com/1657470871/EACKVoLTI?type=comment#_rnd1490771184496
从这开始的内容与标题无关……
之前做过一个小实验：裸写一个单隐层神经网络（实在是很简单），隐层激活函数使用tanh、损失函数使用平方损失函数，样本特征为三维、共三个类，每类10个点。只做训练用，没有测试集。
1. 改变隐层节点个数：使用批处理梯度下降（因为数据量小到可以忽略），迭代500轮后，训练集上的loss变化（图里iters改成epoch更合适）如下图，可以看出来网络表示能力大体上增强（由于随机初始化的原因，每次实验结果可能不同）。
2. 改变学习率：较大的学习率在早期会让loss下降很快，但是后期可能会震荡。
参考资料：
[1] 《神经网络与深度学习讲义》
[2] 《矩阵求导术（上）》
[3] Matrix_calculus