感知机及其代码实现

前言

感知机可以说是一个跨越了人工智能发展历史的模型，其提出之时引起了非常大的关注，但是随后，大家很快发现其存在严重的局限性，这也导致了AI的第一次寒冬。而经过十几年的停滞，多层感知机的出现解决了感知机存在的问题，并且至今仍在被广泛应用，甚至许多前沿的模型架构中都不乏MLP的身影，这也使得感知机成为了学习深度学习一定绕不开的模型，许多经典的深度学习课程与教材都会将感知机作为第一个模型讲解。本文将顺着感知机的发展历史，讲到目前被广泛使用的多层感知机，并介绍与此相关的一些问题，以及相关的一些技术。

注：本文会涉及到梯度下降与线性回归模型的部分内容，对此不了解的读者可以移步本站的梯度下降博客以及线性回归模型博客。

感知机的发展历史

神经元模型

1943年，心理学家McCulloch和数学家Pitts根据生物神经元的结构，提出了MP神经元模型，模型结构如下图

简单来说就是人类的神经元会接受其他神经元的信号，当信号强度超过某一个阈值之后，该神经元被激活。从数学上模拟这个过程就是，多个输入$x_1, x_2, \cdots, x_p$加权相加，这就是收到的信号强度，然后找一个分段函数，$x > 0$给一个激活值，比如定义为$1$，否则不激活，比如定义为$-1$，因为这个函数模拟的是神经元判断信号强度是否达到激活阈值的过程，因此就叫激活函数。整个模型写出来就是

$$y = \phi(w_0 + w_1 x_1 + \cdots + w_p x_p) \\ \phi(x) = \begin{cases} 1 & x \ge 0 \\ -1 & x < 0 \\ \end{cases}$$

由于模型只有一个输出，且激活函数会将输出从连续值映射到$\{-1, 1\}$，因此神经元模型被用于处理二分类问题，这与logistic回归中使用sigmoid函数将连续值映射到$(0, 1)$上的设计是类似的，也可以说，神经元模型等价于换了激活函数的logistic回归。

最后，上面的模型结构图是详细版本，后续的模型随着结构的复杂化，一般都会省略偏置以及损失函数，如下图

感知机模型

在神经元模型的基础上，Frank Roseblatt于1957年提出感知机模型，其模型结构与神经元模型完全相同。

此时大家可能会疑惑，模型结构完全不变，这也能算提出一个新模型？其实不然，感知机的贡献在于其给出了一套机器自主学习参数的算法。生活在机器学习与深度学习盛行年代的我们，可能会先入为主的认为模型自主学习参数是一件很平常的事情，但这种操作最早就是由感知机模型提出的，我们来看看当年感知机模型所设计的算法。

感知机的训练方法为，首先初始化所有参数为$0$，接着取一个样本$(x_1, x_2, \cdots, x_p, y)$，计算得到输出$\hat{y}$，如果$y\hat{y} \le 0$，按下面的公式更新参数

$$w_i \gets w_i - \hat{y} x_i \quad i = 1, 2, \cdots, p \\ w_0 \gets w_0 - \hat{y}$$

否则不更新参数，循环往复。

这个训练方法的想法也很简单，首先，如果预测值是正确的，那么不更新参数，如果预测值是错误的，例如正确值为$1$，预测值为$-1$，此时$w_0 + w_1 x_1 + \cdots + w_p x_p$必然为负，我们希望其最终输出为$1$，那就是希望$w_0 + w_1 x_1 + \cdots + w_p x_p$为正，反映到参数更新上就是我们希望经过参数更新后，$w_0 + w_1 x_1 + \cdots + w_p x_p$的值变大，那只要所有的$w_i x_i$与$w_0$均变大，其和必然变大。综上，我们就得到了参数更新的方向，$w_0$要上调，$w_i$根据$x_i$的正负选择上调或者下调，这样才能保证$w_i x_i$是上调的嘛。现在，再回看上述参数更新公式，就很好理解了。反之亦然，大家可以自行验证。

看完感知机的训练算法，大家可能又会疑惑，模型训练不应该是先选个损失函数，再求解损失函数的最小值吗？而且这个算法怎么看着总有一种梯度下降的既视感？这又是生活在现在的我们先入为主了，在当时是没有把模型训练转化为求解损失函数最小值这种想法的，从算法的设计思路也可以看出来，作者当时就是从结果倒推设计出的这套算法。而看着像梯度下降其实是没错的，从我们现在的角度来看，感知机的这套训练算法等价于定义损失函数为$L(\boldsymbol{w}) = \mathrm{max}(0, -y\hat{y})$，使用随机梯度下降进行优化。

XOR问题

可以进行自主学习的感知机在当时引起了大量关注，不少人因此认为人工智能时代要来了。但是，1969年，Marvin Minsky和Seymour Papert在《Perceptrons》一书中指出，感知机无法解决线性不可分问题，例如XOR问题，即异或问题。XOR问题非常简单，假设有两个输入$x_1$与$x_2$，当两个输入正负相同时则为正类，反之则为负类，等价于对两个输出做异或运算，如下表

$x_1$	$x_2$	$y$
$+$	$+$	$1$
$+$	$-$	$-1$
$-$	$+$	$-1$
$-$	$-$	$1$

一个最简单的XOR数据集只包含四个样本，如下图

其中，不同颜色代表不同类别。现在，我们需要训练一个感知机模型，用于将这四个样本分开。但是，这是做不到的，因为感知机模型本质上还是一个线性模型，因此，这个问题就相当于，找到一条直线，将上图中的四个点分开，很显然这是无法做到的。

XOR问题的解决思路

历史的车轮滚滚向前，在感知机模型被冷落了十几年之后，19世纪80年代，XOR问题也迎来了解决方法。思路很简单，既然一条直线无法将四个点分开，那么就用两个

如上图，如果我们有两个线性模型，那么就可以将XOR问题的四个样本完美的分开。我们来详细看看其决策过程，输入样本$(x_1, x_2)$，如果其位于$f_1$上，或者位于$f_2$下，则该样本属于正类$1$，如果其位于$f_1$下，同时位于$f_2$上，则该样本属于负类$-1$。不难发现，上述决策过程需要将样本同时输入两个线性模型$f_1$与$f_2$，并且整合$f_1$与$f_2$的输出，得到分类结果。因此，我们就得到了下面的模型结构

构建两个线性模型$f_1$与$f_2$，其输出分别为$h_1$与$h_2$，再整合两个模型的输出得到最终输出$y$。从模型结构上来说，就是在感知机的基础上多加了一层中间层，这个中间层就被称为隐含层，而引入了隐含层的感知机就被称为多层感知机。

注：XOR问题只是线性不可分问题中的一个代表，因为其足够简单，只要能够解决XOR问题，那么其他的线性不可分问题也可以被解决。例如，对于多层感知机来说，通过增加隐含层的神经元数量，就可以解决更复杂的线性不可分问题。

多层感知机

感知机的发展历史就先告有一段落，现在我们来详细看看多层感知机（Multi-Layer Perceptron，MLP）模型，模型结构如下图

其中，输入层的节点数由数据的特征数决定，输出层的节点数由具体的问题类型决定，而隐含层的层数与每层的节点数均为超参数，一般均为人为设定。跟上文感知机的模型结构图相同，每一条连接线均代表一个可以学习的参数，且上述结构图中省略了偏置与损失函数。

注1：关于隐含层超参数的设定，初学者常常会困惑：隐含层的神经元数量为什么是$256$？调隐含层的超参数是玄学吗？实际上，隐含层超参数确实没有一套非常具体的调整方法，多以经验为主。一般隐含层的神经元数量会设置为$2$的整数次方，例如$128$，$256$，这主要是因为这样的设置有利于计算机的计算。而隐含层的层数不宜设置过深，过深的MLP在训练中会出现严重的梯度问题。而如果设置了多个隐含层，其神经元数量一般会呈递减，例如输入层神经元数量为$100$，输出层神经元数量为$10$，使用三个隐含层，其神经元数量可以分别设置为$128$，$64$，$32$，尽量避免神经元数量的大幅度变化。而从宏观层面看，更多的隐含层与神经元数量会带来更强的模型复杂度，因此，如果训练结果呈欠拟合，则增加隐含层层数与神经元数量，反之同理。这些基本就是隐含层超参数的调整方法了。

注2：多层感知机还有很多名字，例如前馈神经网络（Feedforward Neural Network，FNN），在神经网络中常常被叫做全连接层，稠密层，线性层，英文中会叫dense layer，大家基本都是按习惯叫法来，在各种资料中看到这些名称的时候知道跟MLP是一个东西就可以了。

注3：感知机模型扩展至多分类问题的方法与softmax回归是一样的OvR策略，将多层感知机的最后一层隐含层到输出层看作是一个softmax回归就很显然了。

激活函数的作用

在感知机模型以及logistic回归中，损失函数的作用是，在处理分类问题中，将连续值映射为离散值或概率值。而在多层感知机中，对于新引入的隐含层，是否需要使用损失函数？为什么要使用损失函数？这是我们这一小节要讲的问题。

先说结论，在多层感知机中，每一层隐含层都必须使用损失函数。至于原因，我们可以反过来想，如果隐含层不使用损失函数会带来什么问题？就用上面所展示的单隐含层的多层感知机，记$\boldsymbol{x} = (x_1, x_2, \cdots, x_p)^T$，$\boldsymbol{h} = (h_1, h_2, \cdots, y_m)^T$，$\boldsymbol{y} = (y_1, y_2, \cdots, y_k)^T$，三层的多层感知机会产生两个参数矩阵，以及两个偏置向量，这里我们选择将权重与偏置分开写，写在一个矩阵里也是没有问题的。从输入层到隐含层的参数矩阵记为$W^{(1)}$，形状为$m \times p$，偏置向量记为$\boldsymbol{w_{0}^{(1)}}$，形状为$m \times 1$，从隐含层到输出层的参数矩阵记为$W^{(2)}$，形状为$k \times m$，偏置向量记为$\boldsymbol{w_{0}^{(2)}}$，形状为$k \times 1$，在不使用损失函数的情况下

$$\boldsymbol{h} = W^{(1)} \boldsymbol{x} + \boldsymbol{w_{0}^{(1)}} \\ \boldsymbol{y} = W^{(2)} \boldsymbol{h} + \boldsymbol{w_{0}^{(2)}}$$

将上式代入下式可以得到

$$\boldsymbol{y} = W^{(2)} (W^{(1)} \boldsymbol{x} + \boldsymbol{w_{0}^{(1)}}) + \boldsymbol{w_{0}^{(2)}} = W^{(2)} W^{(1)} \boldsymbol{x} + (W^{(2)} \boldsymbol{w_{0}^{(1)}} + \boldsymbol{w_{0}^{(2)}})$$

不难发现，即使添加了隐含层，拟合出的模型依旧是一个线性模型，权重为$W^{(2)} W^{(1)}$，偏置为$W^{(2)} \boldsymbol{w_{0}^{(1)}} + \boldsymbol{w_{0}^{(2)}}$，而线性模型是无法解决线性不可分问题的，这与我们引入隐含层的初衷不符。

到这里，我们就可以引出损失函数的另一个作用，为模型增加非线性性。由于单层的感知机是一个线性模型，无法解决线性不可分问题，因此我们选择使用添加隐含层的方法来解决，但是隐含层本身并不具有非线性性，所以隐含层需要使用激活函数来为提供非线性性，隐含层不使用激活函数的多层感知机会退化为单层感知机。

常用激活函数

现在，让我们来总结一下损失函数的作用，第一，在输出层中将连续值映射为离散值或概率值，第二，在隐含层中为模型提供非线性变换。基于此，我们可以得到损失函数需要满足的几个条件

1. 非线性函数
2. 连续可导，但也允许在少数点上不可导
3. 函数尽可能简单
4. 导函数值域在一个合适的区间内，不可过大也不可过小

上述条件中，前两个条件很好理解，用于提供非线性性的函数肯定是非线性的，而连续可导是为了方便模型使用梯度下降进行训练。至于后两个条件，在后文我们会详细讲到，现在，让我们来看几个最常见的损失函数。

sigmoid函数

首先是sigmoid函数，就是logistic回归中使用的激活函数，其表达式为

$$\sigma(x) = \frac{1}{1 + e^{-x}}$$

这也是最经典的损失函数，满足了上述的所有要求，既可以用在隐含层也可以用在输出层，早期的神经网络模型基本均使用sigmoid函数作为损失函数，不过随着技术的发展，现在也逐渐被淘汰了，只会在二分类问题的输出层有所出场。

注：至于softmax函数，基本是为了输出层量身定制的，并不能用在隐含层上，一般做多分类问题的话，输出层就使用softmax函数，这里就顺带提一句了。

tanh函数

然后是tanh函数，这是一个反三角函数，其表达式为

$$\mathrm{tanh}(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}$$

tanh函数可以看作放大并平移的sigmoid函数，由于其值域为$(-1, 1)$，因此一般不会将其用于输出层。提出使用tanh函数来替代sigmoid函数的初衷是为了解决sigmoid函数的导数值域太小的问题，详细原因会在后文关于模型训练的部分细讲。目前，tanh函数的出场空间也很小了，只在RNN中会有一定的使用。

ReLU函数

最后是ReLU（Rectified Linear Unit）函数，也被称为线性整流函数，这是目前使用最广泛的损失函数，没有之一。其表达式为

$$\mathrm{ReLU}(x) = \begin{cases} x & x \ge 0 \\ 0 & x < 0 \\ \end{cases} \\ = \mathrm{max}(0, x)$$

这是一个非常简单的分段函数，而且在$x = 0$处不可导，其值域也不可能被用于输出层上，理应不是一个很好的选择，但是用ReLU函数效果就是好，就这么简单，深度学习很多时候就是这样的。至于效果好的原因，目前普遍认为是因为其$x \ge 0$的部分导数固定为$1$，这对深层的神经网络的训练有很强的正面效果。当然ReLU函数也不是没有缺点，其$x < 0$的部分导数固定为$0$，容易导致部分神经元参数一直得不到更新，即“神经元死亡”，也因此诞生了很多ReLU函数的变体用于缓解该问题，不过后续有dropout技术的出现，ReLU的这个缺点也不是非常明显了，因此目前对于大部分深度神经网络，损失函数的选择依旧是无脑选择ReLU。

注：其实还有大量的基于上述三个损失函数的变体，但是就笔者的个人经历来说，这些变体都属于只在书本上见过，从来没有在实际中用过，因此这里就不讲了，感兴趣的读者可以自行搜索学习。

损失函数与优化算法

最后，就是多层感知机的训练了，到了多层感知机的年代，已经有了将把模型训练转化为求解损失函数最小值的思路了，所以多层感知机的训练与求解线性回归的数值解使用的算法并无二致。

在损失函数的选用上，对于回归问题选用均方误差损失函数（MSE），表达式为

$$L(W) = \sum_{i=1}^{n} (y_i - \hat{y})^2$$

对于分类问题选用交叉熵损失函数，表达式为

$$L(W) = - \sum_{i=1}^{n} \sum_{j=1}^{m} {y}_{i}^{(j)} \mathrm{ln} {\hat{y}}_{i}^{(j)}$$

在优化算法的选择上，均选用梯度下降算法及其变体。

补充：反向传播算法

注：反向传播算法母庸置疑是一个非常重要的算法，但是由于其作用是优化梯度的计算，就是让梯度下降跑得更快，而现在成熟的深度学习框架，如Pytorch，均实现了反向传播算法并将其细节隐藏起来，因此如果是以使用深度学习为目的进行学习，不了解反向传播算法并不影响使用深度学习模型，因此笔者选择将该部分放入补充。但是反向传播算法真的很重要，建议还是看一看。

在学习反向传播算法之前，我们肯定要先知道反向传播算法是干什么的。反向传播算法并不是什么新的优化算法，而是一种用于在深度神经网络中加速梯度计算的算法。上一小节我们就讲到，多层感知机依旧是使用梯度下降及其变体进行优化求解的，那这当中自然会产生大量的梯度计算，而感知机或者说神经网络的层数堆叠，除了带来参数量的剧增，梯度的计算复杂度也会随着增加，尤其是越靠近输入层的参数，其梯度计算越复杂。

而反向传播算法就是说，我们发现参数的梯度计算之间是存在相互依赖的，具体表现为层数靠前的参数的梯度计算会用到层数靠后的参数的梯度计算的中间结果，因此，我们在算梯度的时候，可以从后往前算，然后计算靠后层数的梯度时，将靠前层数需要用到的中间结果储存起来，传给靠前层数使用，这样就可以减少大量的重复计算，由于这种算法需要从后往前计算梯度，因此得名反向传播算法，相对的，计算模型输出的过程则被称之为前向传播。

总结一下就是，原来在梯度下降中算梯度的时候，就是按照每一个参数的梯度计算公式算，但是在深度神经网络中，梯度的计算公式会变得很复杂，计算量会变得很大，而这时候就有人发现，算梯度的过程中是存在很多重复计算的，因此设计了一套算法，然那些被重复计算的结果只需要算一遍，然后存起来复用，这样就大大减少了计算量，这就是反向传播算法。下面，我们会使用一个简单的多层感知机作为例子，从推导每一个参数的梯度计算公式，到发现重复计算的部分，再到最后的反向传播算法。

准备工作

现在，让我们从一个简单的多层感知机模型入手，模型结构与符号标识如下图

从图中可以看到，该模型是一个只有一个隐含层的多层感知机，每一层都只有两个神经元，由于隐含层与输出层都需要经过激活函数，因此使用两个符号分别表示，在隐含层中，用$net_i$表示输入层加权求和的输出结果，用$h_i$表示$net_i$经过激活函数后的结果，例如

$$net_1 = w_{01}^{(1)} + w_{11}^{(1)} x_1 + w_{21}^{(1)} x_2 \\ h_1 = \sigma(net_1)$$

在输出层中同理，用$out_i$与$\hat{y}_i$进行分别表示。从输入层到隐含层的参数记为$w_{ij}^{(1)}$，其中$i$表示输入层的神经元编号，若$i=0$则表示偏置，$j$表示隐含层的神经元编号，例如从$x_1$到$h_1$的参数为$w_{11}^{(1)}$，可以参考上述$net_1$的计算式。同理，从隐含层到输出层的参数记为$w_{ij}^{(2)}$，其中$i$表示隐含层的神经元编号，$j$表示输出层的神经元编号。

最后，为简单起见，我们暂时不考虑小批量的情况，假设只有一个输入样本$((x_1, x_2),(y_1, y_2))$，模型输出值为$(\hat{y}_1, \hat{y}_2)$，隐含层激活函数使用sigmoid函数，输出层激活函数使用softmax函数，损失函数使用交叉熵损失函数。

梯度计算式推导

准备工作结束，接下来我们尝试使用一般的方法来计算参数的梯度，这当中会大量的使用到链式法则。在此之前，我们先把损失函数以及激活函数的导数先算好，以便后续取用

$$L = - \sum_{i=1}^{2} y_i \ln \hat{y}_i, \quad \frac{\partial L}{\partial \hat{y}_i} = -\frac{y_i}{\hat{y}_i} \\ \sigma(x) = \frac{1}{1 + e^{-x}}, \quad \sigma^\prime(x) = \sigma(x)(1 - \sigma(x)) \\ \mathrm{softmax}(x_i) = \frac{e^{x_i}}{\sum\limits_{k} e^{x_k}}, \quad \mathrm{softmax}^\prime(x_i) = \mathrm{softmax}(x_i)(1 - \mathrm{softmax}(x_i))$$

我们先从靠后层数的参数开始计算，以参数$w_{11}^{(2)}$为例，涉及到该参数的前向传播计算有

$$out_1 = w_{01}^{(2)} + w_{11}^{(2)} h_1 + w_{21}^{(2)} h_2 \\ \hat{y}_1 = \mathrm{softmax}(out_1)$$

根据链式法则

$$\frac{\partial L}{\partial w_{11}^{(2)}} = \frac{\partial L}{\partial \hat{y}_1}\frac{\partial \hat{y}_1}{\partial out_1}\frac{\partial out_1}{\partial w_{11}^{(2)}} = -\frac{y_1}{\hat{y}_1} \hat{y}_1(1 - \hat{y}_1) h_1 = -y_1 (1 - \hat{y}_1) h_1$$

同理，其余$w_{ij}^{(2)}$的计算方式是类似的，到这里，其实跟softmax回归中的梯度计算并无区别。接下来，我们来看看靠前层数的参数，以参数$w_{11}^{(1)}$为例，涉及到该参数的前向传播计算有

$$net_1 = w_{01}^{(1)} + w_{11}^{(1)} x_1 + w_{21}^{(1)} x_2, \quad h_1 = \sigma(net_1) \\ out_1 = w_{01}^{(2)} + w_{11}^{(2)} h_1 + w_{21}^{(2)} h_2, \quad \hat{y}_1 = \mathrm{softmax}(out_1) \\ out_2 = w_{02}^{(2)} + w_{12}^{(2)} h_1 + w_{22}^{(2)} h_2, \quad \hat{y}_2 = \mathrm{softmax}(out_2)$$

这里我们可以看到，参数$w_{11}^{(1)}$跟$x_1$与$h_1$有关，而由于层数的堆叠，$h_1$跟后续的所有层的所有中间结果都有关，这也就造成了计算复杂度的剧增，这里的使用的例子还算简单，大家可以尝试一下多加一些层，每一层多加一些神经元，那么除了输出层参数，其余层的参数推导出的梯度计算式都将巨长无比。言归正传，我们还是来看看参数$w_{11}^{(1)}$的梯度计算，同样根据链式法则

$$\frac{\partial L}{\partial w_{11}^{(1)}} = \frac{\partial L}{\partial h_1}\frac{\partial h_1}{\partial net_1}\frac{\partial net_1}{\partial w_{11}^{(1)}} = \frac{\partial L}{\partial h_1} h_1 (1 - h_1) x_1 \\ \frac{\partial L}{\partial h_1} = \frac{\partial L}{\partial \hat{y}_1}\frac{\partial \hat{y}_1}{\partial out_1}\frac{\partial out_1}{\partial h_1} + \frac{\partial L}{\partial \hat{y}_2}\frac{\partial \hat{y}_2}{\partial out_2}\frac{\partial out_2}{\partial h_1} = - y_1 (1 - \hat{y}_1) w_{11}^{(2)} - y_2 (1 - \hat{y}_2) w_{12}^{(2)} \\ \frac{\partial L}{\partial w_{11}^{(1)}} = (y_1 (1 - \hat{y}_1) w_{11}^{(2)} + y_2 (1 - \hat{y}_2) w_{12}^{(2)}) h_1 (h_1 - 1) x_1$$

同理，其余$w_{ij}^{(1)}$的计算方式也是类似的。总结一下就是

$$\frac{\partial L}{\partial w_{ij}^{(2)}} = \frac{\partial L}{\partial \hat{y}_j}\frac{\partial \hat{y}_j}{\partial out_j}\frac{\partial out_j}{\partial w_{ij}^{(2)}} \\ \frac{\partial L}{\partial w_{ij}^{(1)}} = \left( \sum_{k} \frac{\partial L}{\partial \hat{y}_k}\frac{\partial \hat{y}_k}{\partial out_k}\frac{\partial out_k}{\partial h_j} \right) \frac{\partial h_j}{\partial net_j}\frac{\partial net_j}{\partial w_{ij}^{(1)}}$$

梯度计算优化

观察两层参数的计算式，我们不难发现有一些重复计算的部分，例如$\frac{\partial L}{\partial \hat{y}_i}\frac{\partial \hat{y}_i}{\partial out_i}$，那我们不妨记

$$\delta_{i}^{(2)} = \frac{\partial L}{\partial \hat{y}_i}\frac{\partial \hat{y}_i}{\partial out_i}$$

并且不难观察到

$$\frac{\partial out_j}{\partial w_{ij}^{(2)}} = h_i, \quad i \ne 0 \\ \frac{\partial out_j}{\partial w_{ij}^{(2)}} = 1, \quad i = 0$$

故$w_{ij}^{(2)}$的计算式可以写作

$$\frac{\partial L}{\partial w_{ij}^{(2)}} = \delta_{i}^{(2)} h_i, \quad i \ne 0 \\ \frac{\partial L}{\partial w_{ij}^{(2)}} = \delta_{i}^{(2)}, \quad i = 0$$

于是我们可以换一种计算方法，先算出所有的$\delta_{i}^{(2)}$，$h_i$在前向传播中就已经计算完成了，直接拿来用，再根据上述计算式将两者相乘。对比原本的计算方法，即按照计算式以此计算每一个$w_{ij}^{(2)}$的梯度，不难发现，原本的计算方法相当于将每一个$\delta_{i}^{(2)}$都重复计算了$j+1$次，而这还只是在计算所有$w_{ij}^{(2)}$的梯度中节省下来的计算量。

我们继续来看$w_{ij}^{(1)}$的计算式，类比$w_{ij}^{(2)}$的处理方式，观察到

$$\frac{\partial net_j}{\partial w_{ij}^{(1)}} = x_i, \quad i \ne 0 \\ \frac{\partial net_j}{\partial w_{ij}^{(1)}} = 1, \quad i = 0$$

于是，我们也想将$w_{ij}^{(1)}$的计算式转化为如下的形式

$$\frac{\partial L}{\partial w_{ij}^{(1)}} = \delta_{i}^{(1)} x_i, \quad i \ne 0 \\ \frac{\partial L}{\partial w_{ij}^{(1)}} = \delta_{i}^{(1)}, \quad i = 0$$

其中

$$\delta_{i}^{(1)} = \left( \sum_{k} \frac{\partial L}{\partial \hat{y}_k}\frac{\partial \hat{y}_k}{\partial out_k}\frac{\partial out_k}{\partial h_i} \right) \frac{\partial h_i}{\partial net_i}$$

观察发现

$$\frac{\partial L}{\partial \hat{y}_k}\frac{\partial \hat{y}_k}{\partial out_k} = \delta_{k}^{(2)} \\ \frac{\partial out_k}{\partial h_i} = w_{ik}^{(2)}$$

而$\frac{\partial h_i}{\partial net_i}$依损失函数而定，这里就记作

$$\frac{\partial h_i}{\partial net_i} = \phi^\prime(net_i)$$

故

$$\delta_{i}^{(1)} = \left( \sum_{k} \delta_{k}^{(2)} w_{ik}^{(2)} \right) \phi^\prime(net_i)$$

到这里我们就可以发现，正如上文所说，靠前层数的梯度计算是会使用到靠后层数的梯度计算的中间结果的。到了这一层，对比原计算方法，单就$\delta_{k}^{(2)}$的计算次数减少的可就不是一星半点了。

最后，如果再加一层隐含层，梯度计算式会更复杂吗？尝试推导一下就会发现，如果是直接推导计算式，那确实会更复杂，而如果使用上述计算方法，其计算式形式是类似的，假设有一个$\delta_{i}^{(0)}$的话，那其计算式应该为

$$\delta_{i}^{(0)} = \left( \sum_{k} \delta_{k}^{(1)} w_{ik}^{(1)} \right) \phi^\prime(x_i)$$

所以，上面提出的这套算法是可以被应用到任意层数的多层感知机中的，并且$\delta$的计算是靠前层数依赖靠后层数，因此，这套算法的计算需要从后往前算，故得名反向传播算法。

反向传播算法

至此，我们已经通过一个简单的模型了解了反向传播算法是如何优化梯度计算的效率的，现在，我们就给出一般多层感知机模型下的反向传播算法，模型结构图与符号标识如下图

记输出层神经元数量为$p$，隐含层数量为$s$，每一层隐含层的神经元书数量为$m_{i}$，每一层隐含层的两个结果，即经过激活函数前与经过激活函数后的结果，分别用$net^{(i)}_{j}$与$h^{(i)}_{j}$表示，输出层同理，用$out_i$与$\hat{y}_i$进行分别表示，参数则使用$w_{ij}^{(l)}$表示，其中$l$表示参数从前往后的排序，$i$表示上一层的神经元编号，$j$表示下一层的神经元编号，损失函数使用$L$表示，隐含层的激活函数使用$\phi_{i}$表示，输出层的激活函数使用$\phi_{out}$表示。

则使用反向传播算法计算该多层感知机模型的梯度的计算公式为

$$\frac{\partial L}{\partial w_{ij}^{(l)}} = \begin{cases} \delta_{i}^{(l)} h_i & i \ne 0, \ l \ne 1 \\ \delta_{i}^{(l)} x_i & i \ne 0, \ l = 1 \\ \delta_{i}^{(l)} & i = 0 \\ \end{cases} \\ \delta_{i}^{(l)} = \begin{cases} \left( \sum\limits_{k} \delta_{k}^{(l+1)} w_{ik}^{(l+1)} \right) \phi^\prime_{l+1}(net_i^{(l+1)}) & l \ne s+1 \\ \frac{\partial L}{\partial \hat{y}_i}\phi^\prime_{end}(net_i^{(l+1)}) & l = s+1 \\ \end{cases}$$

Pytorch代码实现

现在，我们来介绍一下多层感知机的Pytorch代码实现，运行环境为jupyter notebook

注1：Pytorch是目前最主流的深度学习框架，以下内容将基于此框架实现深度学习算法，对Pytorch不熟悉的读者可以移步Pytorch快速入门作简单了解。

注2：以下内容会介绍如何使用Pytorch框架构建深度学习模型，以及介绍常用的参数，并与上述理论部分相对照，再用Pytorch中自带的小型数据集跑一个简单的demo。

注3：本文并不会将代码讲得面面俱到，在学习Pytorch的过程中建议多多参考PyTorch官方文档或者PyTorch中文文档。

首先，我们需要导入依赖库

import torch
from torch import nn
from torch.utils.data import DataLoader  # 数据读取器
from torchvision import transforms  # 数据预处理
from torchvision.datasets import FashionMNIST  # FashionMNIST数据集

然后，我们需要读取本小节使用的示例数据集，我们将使用FashionMNIST数据集，相信稍微了解深度学习的读者应该都知道MNIST数据集，这是一个经典的手写数字数据集，数据均为灰度图像，类别为10类。但是MNIST数据集对于现代的深度学习模型来说已经过于简单了，因此，FashionMNIST数据集就是为此而生的，其数据为与MNIST数据集形状相同的灰度图像，且也为10个类别的分类问题，可以直接替代MNIST数据集作为图像分类问题的基准数据集。

# 读取数据集
data_FashionMNIST_train = FashionMNIST('.\Data\fashion_mnist_data', train=True, download=True, transform=transforms.ToTensor())
data_FashionMNIST_test = FashionMNIST('.\Data\fashion_mnist_data', train=False, download=True, transform=transforms.ToTensor())

# 生成数据读取器
FashionMNIST_train = DataLoader(data_FashionMNIST_train, batch_size=64, shuffle=True)  # shuffle表示打乱
FashionMNIST_test = DataLoader(data_FashionMNIST_test, batch_size=64, shuffle=True)
# 由于深度学习的数据集一般比较大，因此测试集也会用读取器进行小批量读取，防止爆内存

接着，我们正式介绍多层感知机的实现，首先我们要明白，神经网络模型是由多个网络层堆叠而成，比方说我的第一层是一个线性层，第二层是一个卷积层等等，然后这些层各自会有不同的输入输出神经元数量，常常有人将做神经网络模型比喻为搭积木。Pytorch作为一个深度学习框架，其逻辑就是将常用的网络层一一封装，供用户自行堆叠组装。

前面提到过，多层感知机的每一层在神经网络中通常被称为线性层，那从神经网络的角度来说，多层感知机就是一个由多个线性层堆叠而成的神经网络。Pytorch中用nn.Linear()函数实现线性层，其参数也非常简单，第一个参数为输入神经元数量，第二个参数为输出神经元数量，以FashionMNIST数据集为例，假设我们构建一个两层的多层感知机，由于FashionMNIST数据集中的图像大小为$28 \times 28$，因此第一层的输入为784，而该数据集是一个十个类别的多分类问题，因此第二层的输出为10，再设隐含层大小为256，则第一层的输出与第二层的输入均为256，至此，我们得到了模型的所有参数，构建模型的代码如下

# 多层感知机
class MLP(nn.Module):
    def __init__(self):
        super().__init__()
        self.input = nn.Flatten()
        self.L1 = nn.Linear(784, 256)  # 第一层
        self.L2 = nn.Linear(256, 10)  # 第二层
        self.ReLU = nn.ReLU()  # 激活函数
        
    def forward(self, x):
        x = self.input(x)
        x = self.L1(x)
        x = self.ReLU(x)
        x = self.L2(x)
        return x

将模型实例化，并读取计算设备

MLP = MLP()  # 实例化
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

而后，我们还要解决损失函数，优化算法以及训练，分类问题自然是选用交叉熵损失函数，优化算法选择梯度下降，Pytorch的训练流程相对固定，此处将其封装为一个训练函数，以便重复调用

def model_train(model, train_data, test_data, num_epoch, lr, device):
    # model：模型，train_data：训练数据的DataLoader，test_data：测试数据的DataLoader，num_epoch：epoch数，lr：学习率，device：训练设备
    
    # 定义损失函数与优化算法
    criterion = torch.nn.CrossEntropyLoss()  # 交叉熵
    optimizer = torch.optim.SGD(model.parameters(), lr=lr)  # SGD
    
    model.to(device)  # 将模型放到GPU上
    
    # epoch
    train_loss = []  # 储存每个epoch训练损失
    train_accuracy = []  # 储存每个epoch训练准确度
    test_accuracy = []  # 储存每个epoch测试准确度
    
    for epoch in range(num_epoch):
        # 训练
        model.train()  # 将模型调整为训练模式
        train_run_loss = 0.0  # 储存训练损失
        true_total = 0  # 储存预测正确数量
        data_size = 0  # 储存数据量
        
        for data, target in train_data:
            data, target = data.to(device), target.to(device)  # 将数据放到GPU上
            optimizer.zero_grad()  # 手动梯度归零
            
            output = model(data)  # forward
            loss = criterion(output, target)  # 计算损失值
            loss.backward()  # 反向传播
            optimizer.step()  # 更新参数
            
            train_run_loss += loss.item()  # 累加训练损失，item()方法用于取出loss的数值，损失函数返回的loss不止有数值，还有其他信息
            
            predicted = torch.argmax(output, 1)  # 向量转换为标签
            true_batch = (predicted == target).sum().item()
            true_total += true_batch  # 累加正确数量
            data_size += target.size(0)  # 累加数据量
        
        train_epoch_loss = train_run_loss / len(train_data)  # 本轮epoch的损失值
        train_loss.append(train_epoch_loss)  # 储存本轮epoch损失值
        print('Epoch: {} \tTraining Loss: {:.6f}'.format(epoch+1, train_epoch_loss))
        
        train_epoch_accuracy = true_total / data_size  # 本轮epoch正确率
        train_accuracy.append(train_epoch_accuracy)  # 储存本轮epoch正确率
        print('Epoch: {} \tTraining Accuracy: {:.2f}%'.format(epoch+1, 100 * train_epoch_accuracy))
        
        # 测试
        model.eval()  # 将模型调整为测试模式
        true_total = 0  # 储存预测正确数量
        data_size = 0  # 储存数据量
        with torch.no_grad():  # 关闭梯度计算
            for data, target in test_data:
                data, target = data.to(device), target.to(device)  # 将数据放到GPU上
                
                output = model(data)
                predicted = torch.argmax(output, 1)  # 预测向量转换为标签
                
                true_batch = (predicted == target).sum().item()
                true_total += true_batch  # 累加正确数量
                data_size += target.size(0)  # 累加数据量
                
            test_epoch_accuracy = true_total / data_size  # 本轮epoch正确率
            test_accuracy.append(test_epoch_accuracy)  # 储存本轮epoch正确率
            print('Epoch: {} \tTesting Accuracy: {:.2f}%'.format(epoch+1, 100 * test_epoch_accuracy))
            
    # 绘图
    import matplotlib.pyplot as plt
    plt.figure()
    plt.xlabel("Epoch")
    plt.plot(range(1, num_epoch+1), train_loss, '-', label="Training Loss")
    plt.plot(range(1, num_epoch+1), train_accuracy, '--', label="Training Acc")
    plt.plot(range(1, num_epoch+1), test_accuracy, '-.', label="Testing Accuracy")
    plt.legend()
    
    # 清除GPU缓存
    torch.cuda.empty_cache()

最后，运行训练函数即可

model_train(MLP, FashionMNIST_train, FashionMNIST_test, num_epoch=10, lr=0.1, device=device)

Epoch: 1 	Training Loss: 0.624402
Epoch: 1 	Training Accuracy: 78.22%
Epoch: 1 	Testing Accuracy: 83.09%
Epoch: 2 	Training Loss: 0.437807
Epoch: 2 	Training Accuracy: 84.35%
Epoch: 2 	Testing Accuracy: 83.13%
Epoch: 3 	Training Loss: 0.393465
Epoch: 3 	Training Accuracy: 85.82%
Epoch: 3 	Testing Accuracy: 84.91%
Epoch: 4 	Training Loss: 0.365517
Epoch: 4 	Training Accuracy: 86.64%
Epoch: 4 	Testing Accuracy: 86.37%
Epoch: 5 	Training Loss: 0.345038
Epoch: 5 	Training Accuracy: 87.56%
Epoch: 5 	Testing Accuracy: 86.41%
Epoch: 6 	Training Loss: 0.330205
Epoch: 6 	Training Accuracy: 88.05%
Epoch: 6 	Testing Accuracy: 85.13%
Epoch: 7 	Training Loss: 0.316345
Epoch: 7 	Training Accuracy: 88.55%
Epoch: 7 	Testing Accuracy: 83.08%
Epoch: 8 	Training Loss: 0.303760
Epoch: 8 	Training Accuracy: 88.86%
Epoch: 8 	Testing Accuracy: 85.90%
Epoch: 9 	Training Loss: 0.296070
Epoch: 9 	Training Accuracy: 89.17%
Epoch: 9 	Testing Accuracy: 86.96%
Epoch: 10 	Training Loss: 0.286688
Epoch: 10 	Training Accuracy: 89.58%
Epoch: 10 	Testing Accuracy: 86.72%

正则化与丢弃法

了解完多层感知机后，现在我们来介绍一下在多层感知机的训练中可能会遇到的问题，以及相关的一些技术，例如本小节我们将讨论的过拟合问题，以及相关的技术正则化与丢弃法。

正则化

对于正则化，相信接触过线性回归的读者应该并不陌生，其通过在损失函数中添加正则化项，以限制模型的复杂度，从而实现防止过拟合的目的。常用的正则化项有$L_1$正则化与$L_2$正则化，$L_1$正则化项就是模型的所有参数的绝对值之和，而$L_2$正则化就是模型的所有参数的平方和，最后，正则化项还会乘以一个正则化系数$\alpha$，该系数用于控制正则化的强度，系数越大，正则化强度越大，模型的复杂度越低，反之同理。而这一套方法在多层感知机模型中也是同样适用的，在损失函数后添加正则化项，而后正常进行训练即可，例如MSE损失函数，添加$L_2$正则化后如下

$$L(W) = - \sum_{i=1}^{n} \sum_{j=1}^{m} {y}_{i}^{(j)} \mathrm{ln} {\hat{y}}_{i}^{(j)} + \alpha \Omega(W)$$

其中，$\Omega(W)$为模型的所有参数的平方和，即正则化项。

丢弃法

而丢弃法（Dropout）是一种针对多层神经网络结构的模型设计的，用于防止过拟合问题的方法，对于神经网络模型来说，相比于正则化，Dropout是更常用的防止过拟合的方法。

Dropout的策略是对于隐含层（非输入输出层），在训练过程中，对每一轮迭代（即每一个小批量的训练），均按相同的概率$p$随机丢弃部分神经元，再进行训练，如下图

其实现方法为，生成一个与隐含层输出向量$\boldsymbol{h}$形状相同的随机掩码（Mask）向量$\boldsymbol{m}$，该向量中的每一个元素有$p$的概率为$0$，$1 - p$的概率为$1$，将隐含层的输出向量$\boldsymbol{h}$与掩码（Mask）向量$\boldsymbol{m}$按元素相乘，则被丢弃的神经元的输出会变为$0$，此时该神经元对后续的网络层均不会产生任何影响，也就相当于该神经元被丢弃。

在模型的推理阶段，Dropout将不再丢弃神经元，但是此时会出现一个问题，由于在训练过程中仅有部分神经元参与工作，如果在推理过程中启用所有的神经元，那么可以预见在训练过程中使用了Dropout的网络层的输出会变大，这显然会出大问题，但是如果在推理过程中也随机进行神经元的丢弃，那么会给模型引入随机性，这也是我们不愿看到的。

因此，最终的解决方案为，在推理阶段，对使用了Dropout的网络层的输出进行缩放，具体为将每一个输出值乘以$1 - p$，其实不难理解，假设一个隐含层的输出之和的期望值为$x$，在训练阶段其仅启用$1 - p$的神经元，因此训练阶段的输出期望为$x(1 - p)$，而推理阶段会启用全部的神经元，此时输出期望为$x$，要保证训练阶段与推理阶段的输出期望一致，只需将推理阶段的输出值乘以$1 - p$即可。

Pytorch代码实现

对于上述介绍的技术，Pytorch也提供了相应的实现，首先是正则化，Pytorch将$L_2$正则化的实现封装在优化算法中，在所有的优化算法中均有一个参数weight_decay，该参数即为正则化系数$\alpha$，默认为$0$，即不使用正则化。而$L_1$正则化在Pytorch中并没有被实现，若需使用，则必须进行手动实现。在上述训练函数中，将优化算法对应的代码改为如下代码，并为训练函数添加一个新的参数wd，即可在模型训练中使用$L_2$正则化

torch.optim.SGD(model.parameters(), lr=lr, weight_decay=wd)

而Dropout在Pytorch中被视为一个网络层，其使用nn.Dropout()函数实现，该函数有一个参数p，对应Dropout中的丢弃概率$p$。该网络层的作用为有概率$p$将输入变为$0$，因此，只需要将其放在使用Dropout的网络层后即可，例如，对上文中实现的多层感知机模型的隐含层使用Dropout，代码如下

# 丢弃法（等价于正则化）
class MLP_dropout(nn.Module):
    def __init__(self):
        super().__init__()
        self.input = nn.Flatten()
        self.L1 = nn.Linear(784, 256)
        self.L2 = nn.Linear(256, 10)
        self.ReLU = nn.ReLU()
        self.Dropout = nn.Dropout(p=0.5)
        
    def forward(self, x):
        x = self.input(x)
        x = self.L1(x)
        x = self.ReLU(x)
        x = self.Dropout(x)
        x = self.L2(x)
        return x

MLP_dropout = MLP_dropout()
model_train(MLP_dropout, FashionMNIST_train, FashionMNIST_test, num_epoch=10, lr=0.1, wd=0, device=device)

Epoch: 1 	Training Loss: 0.681472
Epoch: 1 	Training Accuracy: 76.13%
Epoch: 1 	Testing Accuracy: 82.08%
Epoch: 2 	Training Loss: 0.484318
Epoch: 2 	Training Accuracy: 82.67%
Epoch: 2 	Testing Accuracy: 81.58%
Epoch: 3 	Training Loss: 0.441137
Epoch: 3 	Training Accuracy: 84.11%
Epoch: 3 	Testing Accuracy: 83.53%
Epoch: 4 	Training Loss: 0.415889
Epoch: 4 	Training Accuracy: 85.00%
Epoch: 4 	Testing Accuracy: 85.39%
Epoch: 5 	Training Loss: 0.397788
Epoch: 5 	Training Accuracy: 85.68%
Epoch: 5 	Testing Accuracy: 84.75%
Epoch: 6 	Training Loss: 0.386381
Epoch: 6 	Training Accuracy: 86.04%
Epoch: 6 	Testing Accuracy: 86.20%
Epoch: 7 	Training Loss: 0.374492
Epoch: 7 	Training Accuracy: 86.47%
Epoch: 7 	Testing Accuracy: 86.19%
Epoch: 8 	Training Loss: 0.366260
Epoch: 8 	Training Accuracy: 86.62%
Epoch: 8 	Testing Accuracy: 86.85%
Epoch: 9 	Training Loss: 0.356786
Epoch: 9 	Training Accuracy: 86.93%
Epoch: 9 	Testing Accuracy: 87.07%
Epoch: 10 	Training Loss: 0.350871
Epoch: 10 	Training Accuracy: 87.34%
Epoch: 10 	Testing Accuracy: 85.82%

从结果上我们可以看出，相比于不使用正则化与丢弃法的结果，训练准确率与测试准确率之间的差距明显更小。

参数初始化

最后，我们来讨论深度学习模型训练中的另一个重要问题，参数初始化问题，顾名思义，就是训练前，深度学习模型的初始权重应该如何设置。千万不要觉得参数初始化就是随便生成一些随机数，一个好的初始参数可以让模型的训练事半功倍，甚至可能影响模型的最终效果。

对称性问题

一般来说，我们会希望初始参数在$0$附近，这是因为激活函数一般在$0$附近会有较大的梯度，这更有利于模型训练初始阶段的快速收敛。而本小节，我们先来回答一个初学者很容易有疑问的问题，为什么不能初始化参数为$0$？这涉及到著名的对称性问题，这里我们沿用上文中反向传播的例子与符号

我们关注$w_{11}^{(2)}$与$w_{21}^{(2)}$两个参数，这两个参数的梯度计算式为

$$\frac{\partial L}{\partial w_{11}^{(2)}} = -y_1 (1 - \hat{y}_1) h_1 \\ \frac{\partial L}{\partial w_{21}^{(2)}} = -y_1 (1 - \hat{y}_1) h_2$$

由于所有的参数均被初始化为$0$，因此第一轮训练时，前向传播计算出的$h_1$与$h_2$是相等的，而根据上述梯度计算式，$w_{11}^{(2)}$，$w_{21}^{(2)}$两个参数会得到相同的更新，其他参数同理，这又会导致下一轮训练的前向传播计算出的$h_1$与$h_2$也相等，如此往复，这就是著名的对称性问题，显然，这样训练出的模型肯定是有问题的。

梯度爆炸与梯度消失

让我们回到正题，我们希望初始参数在$0$的附近，那么很自然可以想到选择一个均值为$0$，方差很小的概率分布，生成服从该分布的随机数，作为模型的初始参数，例如$N(0, 0.01)$，对于简单的神经网络模型这是可行的。

由于现代的神经网络模型，其网络深度都很深，从上文我们可以得知，梯度下降的过程中，浅层梯度的计算式中存在深层梯度的连乘，在网络层数足够深时，这种计算方法就很容易造成浅层梯度的数值过大或者过小，这也就是著名的梯度爆炸与梯度消失问题。

解决梯度爆炸与梯度消失的方法有很多，例如梯度裁剪，通过网络结构设计在乘法运算中加入加法运算（ResNet），而合适的参数初始化也是解决方法之一。

方差守恒与Xavier初始化

刚才说到，我们会选择一个均值为$0$，方差很小的概率分布用于进行参数初始化，均值已经确定，那我们自然要在方差上做文章。这里我们就要引入一个重要的支撑理论——方差守恒，其想法为，若在每一个网络层中，在正向传播与反向传播的过程中，其输入与输出值的方差相等，则可以避免训练中出现梯度爆炸与梯度消失。

这样说可能有点抽象，下面我们通过Xavier初始化的推导来详细了解如何实现方差守恒，Xavier初始化是一种经典且实用的初始化，从经验上讲，Xavier初始化可以应对大多数情况，也是大多数教材会介绍的第一种初始化方法。

下面我们正式开始推导，首先，我们要聚焦于一个网络层，不考虑激活函数，记其输入向量的元素为$x_j$，向量维度为$n_{in}$，输出向量元素为$h_i$，向量维度为$n_{out}$，权重记为$w_{ij}$，假设$w_{ij}$均值为$0$，方差为$\sigma^2$，$x_j$均值为$0$，方差为$\mathrm{Var}(x_j)$，$w_{ij}$与$x_j$均独立同分布且相互独立，该网络层的输入$x_j$与输出$h_i$之间的计算式为

$$h_i = \sum_{j = 1}^{n_{in}} w_{ij} x_j$$

则$h_i$均值也为$0$，方差记为$\mathrm{Var}(h_i)$，根据方差守恒，我们希望输入与输出的方差相等，即$\mathrm{Var}(h_i) = \mathrm{Var}(x_j)$，我们尝试计算$\mathrm{Var}(h_i)$

$$\mathrm{Var}(h_i) = \mathrm{E}(h_i^2) - [\mathrm{E}(h_i)]^2 = \mathrm{E}(h_i^2) = \mathrm{E}\left[ (\sum_{j = 1}^{n_{in}} w_{ij} x_j)^2 \right] = \mathrm{E}\left[ \sum_{j = 1}^{n_{in}} w_{ij}^2 x_j^2 + \sum_{j \neq k} w_{ij} w_{ik} x_j x_k \right] \\ = \mathrm{E}\left[ \sum_{j = 1}^{n_{in}} w_{ij}^2 x_j^2 \right] = \sum_{j = 1}^{n_{in}} \mathrm{E}(w_{ij}^2) \mathrm{E}(x_j^2) = \sum_{j = 1}^{n_{in}} \mathrm{Var}(w_{ij}) \mathrm{Var}(x_j) = n_{in} \sigma^2 \mathrm{Var}(x_j)$$

注：上述推导中多处使用了$0$均值以及独立性进行计算，且多处使用了公式$\mathrm{Var}(x) = \mathrm{E}(x) - [\mathrm{E}(x)]^2$

经过上述推导，我们可以得到结论，若要满足目标$\mathrm{Var}(h_i) = \mathrm{Var}(x_j)$，则需要满足$n_{in} \sigma^2 = 1$。

不过还没完，方差守恒理论中说的是正向传播与反向传播都要满足方差相等，因此我们还需要计算反向传播的情况，依旧沿用上述符号以及假设，在反向传播中，输入与输出跟正向传播是相反的，且均为梯度值，即输入为${\partial L}/{\partial h_i}$，输出为${\partial L}/{\partial x_j}$，两者计算式为

$$\frac{\partial L}{\partial x_j} = \sum_{i = 1}^{n_{out}} w_{ij} \frac{\partial L}{\partial h_i}$$

我们希望$\mathrm{Var}({\partial L}/{\partial x_j}) = \mathrm{Var}({\partial L}/{\partial h_i})$，同样尝试计算$\mathrm{Var}({\partial L}/{\partial x_j})$

$$\mathrm{Var} \left( \frac{\partial L}{\partial x_j} \right) = \mathrm{E} \left[ \left( \frac{\partial L}{\partial x_j} \right)^2 \right] - \left[ \mathrm{E}\left( \frac{\partial L}{\partial x_j} \right) \right]^2 = \mathrm{E} \left[ \left( \sum_{i = 1}^{n_{out}} w_{ij} \frac{\partial L}{\partial h_i} \right)^2 \right] - \left[ \mathrm{E}\left( \sum_{i = 1}^{n_{out}} w_{ij} \frac{\partial L}{\partial h_i} \right) \right]^2 \\ = \mathrm{E} \left[ \left( \sum_{i = 1}^{n_{out}} w_{ij} \frac{\partial L}{\partial h_i} \right)^2 \right] = \sum_{i = 1}^{n_{out}} \mathrm{E}(w_{ij}^2) \mathrm{E} \left[ \left( \frac{\partial L}{\partial h_i} \right)^2 \right] \\ = \sum_{i = 1}^{n_{out}} \mathrm{Var}(w_{ij}) \mathrm{Var} \left( \frac{\partial L}{\partial h_i} \right) = n_{out} \sigma^2 \mathrm{Var} \left( \frac{\partial L}{\partial h_i} \right)$$

计算过程与正向传播非常相似，而得出的结论为，若要满足$\mathrm{Var}({\partial L}/{\partial x_j}) = \mathrm{Var}({\partial L}/{\partial h_i})$，则需要满足$n_{out} \sigma^2 = 1$，显然，两个条件无法同时被满足，因此，我们会选择满足一个折中的条件$\frac{1}{2}(n_{in} + n_{out}) \sigma^2 = 1$。

至此，我们就完成了对Xavier初始化的推导，先别着急问这个结论有什么用，让我们先来回顾一下我们干了什么。首先，我们在寻找一个参数初始化的方法，也就是给上述推导中的$w_{ij}$找一个初始值，然后，我们想到了选择一个均值为$0$，方差未定的概率分布，从中抽值作为初始值，那剩下的问题就是如何确定方差，也就是上述推导中的$\sigma^2$，而在上述推导中，可以得到

$$\frac{1}{2}(n_{in} + n_{out}) \sigma^2 = 1 \quad \Rightarrow \quad \sigma^2 = \frac{2}{n_{in} + n_{out}}$$

也就是说，Xavier初始化就是帮我们确定了这个未定的方差值，需要注意的是，这个方差值与网络层的输入输出数量有关，也就是说，在Xavier初始化中，不同网络层所使用的概率分布的方差是不同的，以上文实现的多层感知机为例，一共有两个网络层，第一个网络层中$n_{in} = 784$，$n_{out} = 256$，第二层中$n_{in} = 256$，$n_{out} = 10$，两层进行参数初始化时，采用的分布是不同的。

最后，还剩最后一件事情，就是确定分布类型，这里一般使用正态分布或者均匀分布，Pytorch中的实现会使用均匀分布，具体原因有很多，比如正态分布有概率出现极端数值，均匀分布产生的随机数更便于计算等等，不过从实际使用上来讲，两者差别其实不大，根据Xavier初始化的均值和方差，使用正态分布与均匀分布进行初始化时所使用的具体分布为

$$N \left(0, \sqrt{\frac{2}{n_{in} + n_{out}}} \right) \quad or \quad U \left( -\sqrt{\frac{6}{n_{in} + n_{out}}}, \sqrt{\frac{6}{n_{in} + n_{out}}} \right)$$

Pytorch代码实现

Pytorch自然也实现了多数常用的参数初始化方法，下面我们就来介绍一下，在Pytorch中实现参数初始化一般会采用自定义初始化函数配合模型的apply()方法使用，还是以上文中实现的多层感知机为例，参数初始化函数如下

# 参数初始化
def weight_init(m):  # 初始化函数，对不同的层使用不同的初始化方法
    # 线性层参数初始化
	if isinstance(m, nn.Linear):
    	nn.init.xavier_normal_(m.weight)
        nn.init.constant_(m.bias, 0)

该函数中，输入m为神经网络模型的网络层，isinstance()函数用于判断该网络层是否为指定类型，例如此处判断网络层是否为线性层。常用的参数初始化方法均在nn.init中，此处用到了nn.init.xavier_normal_()，即Xavier初始化，输入为网络层的权重，即代表对网络层的权重使用Xavier初始化，还用到了nn.init.constant_()，该函数用于将参数初始化为常数，此处代表将网络层的偏置初始化为$0$。

而后，我们使用模型的apply()方法，参数为上述初始化函数名

model.apply(weight_init)  # 使用model的apply方法，将初始化函数递归应用于网络的所有层

在Pytorch中，模型的apply()方法的作用是，遍历模型的所有网络层，并将指定函数应用于模型的网络层，上述代码的作用就是将实现好的参数初始化函数遍历的作用于模型的所有网络层。

将上述参数初始化代码加入到训练函数中，再进行一次训练

model_train(MLP, FashionMNIST_train, FashionMNIST_test, num_epoch=10, lr=0.1, wd=0, device=device)

Epoch: 1 	Training Loss: 0.573314
Epoch: 1 	Training Accuracy: 79.96%
Epoch: 1 	Testing Accuracy: 79.53%
Epoch: 2 	Training Loss: 0.419421
Epoch: 2 	Training Accuracy: 85.01%
Epoch: 2 	Testing Accuracy: 84.63%
Epoch: 3 	Training Loss: 0.376108
Epoch: 3 	Training Accuracy: 86.46%
Epoch: 3 	Testing Accuracy: 85.10%
Epoch: 4 	Training Loss: 0.350751
Epoch: 4 	Training Accuracy: 87.40%
Epoch: 4 	Testing Accuracy: 86.14%
Epoch: 5 	Training Loss: 0.333174
Epoch: 5 	Training Accuracy: 87.85%
Epoch: 5 	Testing Accuracy: 86.60%
Epoch: 6 	Training Loss: 0.317027
Epoch: 6 	Training Accuracy: 88.35%
Epoch: 6 	Testing Accuracy: 87.00%
Epoch: 7 	Training Loss: 0.304352
Epoch: 7 	Training Accuracy: 88.90%
Epoch: 7 	Testing Accuracy: 87.39%
Epoch: 8 	Training Loss: 0.293761
Epoch: 8 	Training Accuracy: 89.31%
Epoch: 8 	Testing Accuracy: 85.86%
Epoch: 9 	Training Loss: 0.284190
Epoch: 9 	Training Accuracy: 89.53%
Epoch: 9 	Testing Accuracy: 87.42%
Epoch: 10 	Training Loss: 0.277256
Epoch: 10 	Training Accuracy: 89.74%
Epoch: 10 	Testing Accuracy: 87.37%

可以看出，相比与不使用Xavier初始化，第一轮训练的损失值更小，且后续模型的收敛速度也更快。