Book Notes on Pattern Recognition and Machine Learning (PRML)
在之前的线性模型中，我们的模型基于固定非线性基函数 \(\phi(x)\) 的线性组合，即
\[y(x) = f(\sum_{i=1}^M w_i\phi_i(x))\]在分类任务中，\(f(\cdot)\) 为非线性激活函数；在回归任务中，\(f(x)=x\) 。接下来我们进一步考虑，如果不用固定的基函数 \(\phi(x)\) ，而是让基函数内的参数也可以修改调整，这就引出了神经网络。
在最基础的神经网络中，每一个基函数本身又可以看作是输入的线性组合。我们首先考虑最简单的两层神经网络。我们假设输入为 \(\{x_0,...,x_D\}\) ，中间层为 \(\{h_0,...,h_M\}\) ，最后的输出为 \(\{y_1,...,y_K\}\) （零下标的引入表示偏置bias）
在第一层神经网络中，网络的输入为输入数据 X，输出为中间层数据 H，此时可以写出这一层对应的关系
\[h_m = f(\sum_{i=0}^D w_{mi}^{(1)}x_i)\]在第二层神经网络中，网络的输入为中间层数据 H，输出为最终结果 Y，同样可以写出对应关系
\[y_k = g(\sum_{i=0}^Mw_{ki}^{(2)}h_i)\]我们把两层网络结合到一起，即可得到数据输出与输入之间的关系
\[y_k(x) = g(\sum_{i=0}^Mw_{ki}^{(2)}f(\sum_{j=0}^D w_{ij}^{(1)}x_j))\]这里的 \(f(\cdot)\) 和 \(g(\cdot)\) 都代表激活函数。通常选用连续的 sigmoidal 函数，利于整个网络微分，具体的激活函数选择视具体问题确定。但不管怎样，激活函数都必须具有非线性的性质，否则多层的网络就失去了意义，多重线性叠加等价于单层线性。
在实际中，可以用更多层叠加构成网络，此时包括一个输入层，一个输出层和多个中间层
\[h_l=f_l(s_l), \quad s_l = w_lh_{l-1}+b_l,\quad for\ l=1,2,...,L\]这里 \(l\) 代表不同的层，\(h_0\) 代表输入 X，\(h_L\) 用于预测输出 Y
需要注意的是最后一层输出后，还需要经过一个激活函数得到最后的预测结果。这个激活函数记为输出激活函数。在不同任务中，使用不同的输出激活函数并对应不同的loss。在回归任务中，我们不需要输出激活函数，直接使用最小二乘loss即可。在分类任务中，我们最后要给出的是各个类别的分类预测概率，所以需要输出激活函数需要用 logistics sigmoid(二分类)或 softmax(多分类)，并结合交叉熵Loss进行优化。
现在我们已经完成了多层感知机的建模，下一步就是考虑如何去训练优化我们的模型。我们用 \(w\) 统一代表模型所有的待优化参数，例如两层感知机中的 \(w^{(1)}, w^{(2)}\) 。我们的目的是找到一个 \(w\) 使得对应的 Loss 值最小。在这个问题中，很难给出梯度为零的模型解析解，因此通常采用迭代更新的方法逐渐修改模型参数，即
\[w^{(t+1)} = w^{(t)} + \Delta w^{(t)}\]而在迭代更新中，我们可以结合梯度信息加速迭代到损失函数极小值的过程，使得每次参数都沿着梯度下降的方向更新，即梯度下降法。
我们用 \(L\) 来表示 Loss function 损失函数。最直接的利用梯度信息更新参数的方法就是
\[w^{(t+1)} = w^{(t)} - \eta \nabla L(w^{(t)})\]这里的 \(\eta\) 为正数，表示学习率，控制更新的幅度。但这种方法的更新速度显然太慢了，每次更新参数都需要计算整个数据集的 Loss 再求梯度。为此出现了一种在线算法，称为随机梯度下降，每次只选择数据集中的一个数据求其 Loss 的梯度用来更新，即
\[w^{(t+1)} = w^{(t)} - \eta \nabla L_i(w^{(t)}), \quad where\ \sum_{i=1}^n L_i = L\]这里比较一下两种方法的优劣
为了平衡梯度与效率，我们可以选择折中的方式，将数据集划分为一个个batch，每次随机选择一个batch，用该batch内的数据计算梯度，更新模型参数。
在我们引入了模型参数迭代更新方法后，下一步就是求损失函数的梯度信息。在神经网络中，梯度信息是通过误差反向传播来高效求解的。
这里我们继续沿用上面的两层网络的例子。为了区分训练过程中网络的预测输出值和ground-truth标签值，我们用 P 表示网络的预测，用 Y 表示 ground-truth，其他定义和上面一样，这里复制一下。假设输入为 \(\{x_0,...,x_D\}\) ，中间层为 \(\{h_0,...,h_M\}\) （零下标的引入表示偏置bias）
\[p_k(x) = g(\sum_{i=0}^Mw_{ki}^{(2)}h_i) = g(\sum_{i=0}^Mw_{ki}^{(2)}f(\sum_{j=0}^D w_{ij}^{(1)}x_j))\]为了简化分析，我们假设激活函数均采用 logistics sigmoid ，并且是二分类任务。我们可以写出极大似然的目标
\[max \log P = \log \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} = \sum_{i=1}^n [y_iw^{(2)}h-\log (1+e^{w^{(2)}h})]\]将其展开并转换为 Loss 可以得到
\[L(w^{(1)},w^{(2)}) = \sum_{i=1}^n \{y_i \sum_{j=1}^M w^{(2)}_{ij}h_j-\log (1+exp(\sum_{j=1}^M w^{(2)}_{ij}h_j)\}\]我们首先对 \(w^{(2)}\) 求导可得
\[\frac{\partial L}{\partial w^{(2)}_{ij}} = (y_i-p_i)h_j\]进一步根据链式法则对 \(w^{(1)}\) 求导可得
\[\frac{\partial L}{\partial w^{(1)}_{jk}} = \frac{\partial L}{\partial h_j}\frac{\partial h_j}{\partial w^{(1)}_{jk}}\]L 对 h 的求导与对 \(w^{(2)}\) 的求导类似，我们直接写出
\[\frac{\partial L}{\partial h_j} =\sum_{i=1}^n (y_i-p_i)w^{(2)}_{ij}\]而 h 对 \(w^{(1)}\) 的求导，我们根据 sigmoid 求导可得
\[\frac{\partial h_j}{\partial w^{(1)}_{jk}} = -\frac{exp(-w^{(1)}x_j) (-x_j)}{(1+exp(-w^{(1)}x_j))^2} = h_j(1-h_j)x_k\]综合以上两个求导，我们得到最终 L 对 \(w^{(1)}\) 的求导结果
\[\frac{\partial L}{\partial w^{(1)}_{jk}} = \sum_{i=1}^n (y_i-p_i)w^{(2)}_{ij}h_j(1-h_j)x_k\]由此我们也可以看出反向传播这个命名的来源。相比于从输入数据到预测输出的正向传播，我们在求梯度时，信息从输出层经过中间层传至输入层。我们首先根据预测输出和标签值，得到最后一层网络的梯度；再以此类推，直到得到第一层网络的梯度。
在反向传播的第一步，我们需要传递的是 \(y-p\) 即误差。这就要求我们的输出层激活函数必须是 sigmoid 或 softmax，只有这样我们才能得到导数表达式中的 error 项（分类任务）。而中间层的激活函数没有要求，就只有求导计算式的不同，还是以上面的例子说明
\[Tanh: f(x) = \frac{e^x-e^{-x}}{e^x+e^{-x}} \Rightarrow \frac{\partial h_j}{\partial w^{(1)}_{jk}} = (1-h_j^2) x_k\] \[ReLU: f(x) = max(0,x) \Rightarrow \frac{\partial h_j}{\partial w^{(1)}_{jk}} = 1(h_j>0)\cdot x_k\]最后我们写一下反向传播的矩阵形式。首先写出链式法则的矩阵形式。假设 \(Y = (y_i)_{m\times 1}\) ，\(X = (x_i)_{n\times 1}\)
如果 \(Y=g(X)\) ，则可定义
如果 \(Y=g(H), H=f(X)\) ，则有
\[\frac{\partial y_i}{\partial x_j} = \sum_{m} \frac{\partial y_i}{\partial h_m}\frac{\partial h_m}{\partial x_j} \Rightarrow \frac{\partial Y}{\partial X^T} = \frac{\partial Y}{\partial H^T}\frac{\partial H}{\partial X^T}\]这里我们保持和之前多层感知机的例子一致，即
\[h_l=f_l(s_l), \quad s_l = w_lh_{l-1}+b_l,\quad for\ l=1,2,...,L\]这里 \(l\) 代表不同的层，\(h_0\) 代表输入 X，\(h_L\) 用于预测输出 Y
我们首先写出链式法则的关系
我们可以得到最后的结论
\[\frac{\partial L}{\partial h_{l-1}} = W_l^Tf_l'\frac{\partial L}{\partial h_l}, \qquad \frac{\partial L}{\partial W_l} = f_l'\frac{\partial L}{\partial h_l}h_{l-1}^T\]这里的 \(f_l' = \frac{\partial h_l}{\partial s_l^T}\) 是对角矩阵，第 k 个元素是 \(\frac{\partial h_{lk}}{\partial s_{lk}}\)
我们首先给出 Hessian Matrix 的定义。我们将 Loss function L(w) 在某一点 \(\hat{w}\) 泰勒展开可以得到
\[L(w) \approx L(\hat{w}) + (w-\hat{w})^T\nabla L|_{w=\hat{w}} + \frac{1}{2} (w-\hat{w})^T H (w-\hat{w})\]这里的 H 即为 Hessian Matrix，\(H = \nabla \nabla L\) ，展开可以写成
\[H_{ij} = \frac{\partial^2 L}{\partial w_i \partial w_j}\vert_{w=\hat{w}}\]接下里我们考虑用反向传播的思想来求 Hessian Matrix ，PRML书中提到了几种求 Hessian Matrix 的近似方法
这种方法只适用于训练好的网络，因为我们在近似过程中忽略的一项与误差相关。在更广泛的网络下，误差会很大，此时这一项则不能忽略，这种近似变得不精确
所以我们假设初始矩阵 \(H_0 = \alpha I\)，这里 \(\alpha\) 是一个小量，最后就可以根据这个迭代算法找到 \(H+\alpha I\) 的逆
Hessian 矩阵的作用有很多，接下来举一个例子来说明其重要性。
在我们训练网络的时候，如果某一时刻 Loss 的梯度突然为零。此时按照 SGD 的迭代方式，模型参数将不再更新。此时有两种情况，一种是我们来到了局部的极小值，另一种是我们来到了鞍点 saddle point
那么如何区分具体是哪一种情况呢，我们需要利用 Hessian Matrix 来了解 Loss surface 的情况。这里直接给出结论
如果来到了鞍点，我们可以也利用 Hessian Matrix 来进一步确定参数更新的方向，从而逃离鞍点。我们取小于零的 eigen value 值对应的 eigen vector 的方向作为参数更新方向。（这种方法在实际中很少用，因为涉及 H 的计算量很大）
随着维度的增加，在低维中看起来是局部极小值的点，可能在高维中仍然可以优化。所以当模型参数很多，形成的维度空间很大时，local minima 可能本身就是很罕见的。
在神经网络模型中，我们需要控制模型复杂度来避免过拟合的现象。控制模型复杂度有很多方法，其中一种方法是选择一个较大的模型参数数量，在损失函数中添加正则项来控制复杂度。最简单的正则项是二次函数，也被称为权值衰减，即
\[\tilde{L}(w) = L(w) + \lambda w^Tw\]这种正则方式在之前的线性模型中介绍过，超参 \(\lambda\) 控制模型的复杂度。但这种正则方式不能在线性缩放中保持一致。以两层感知机为例，如果我们同时以不同倍数缩放输入值和标签值，那么两层网络参数的变化倍数是不同的，而仅凭 \(\lambda w^Tw\) 这一项则无法进行区分。
此外，另一种控制模型复杂度的方法称为早停止 early stopping，是说我们在训练的过程中，添加一个validation set，在validation set上同样去计算误差，当误差达到最小时就提前停止训练。
在模式识别中，我们的模型需要具有 Invariances 的性质。也就是说，当我们对输入数据进行变换后，对应的输出预测结果应该是不变的。首先我们用数学来刻画表征这样的变换。
以图像数据为例，假设 x 是一张输入图片，s 表示变换（比如旋转），并由一个变换系数 \(\xi\) 控制（比如旋转角度）。此时变换后的数据可表示为 \(s(x,\xi)\) ，同时 \(s(x,0) = x\) 表示无变换。
接下来我们介绍两种增加模型泛化能力的方法。第一种方式是最直觉的，我们直接对输入数据进行不同程度的随机变换，将变换后的数据输入模型进行训练。以最小二乘误差为例，原本的 Loss 可写为
\[L = \int\int (y(x)-y^*)^2p(y^*|x)p(x)dxdy^*\]而对其变换后可写为
\[L = \int\int \{y(s(x,\xi))-y^*\}^2p(y^*|x)p(x)p(\xi)dxdy^*d\xi\]第二种方式是通过添加正则项的方式，这种方法也称为 tangent propagation，我们计算出变换函数的导数 \(\tau\)
\[\tau = \frac{\partial s(x,\xi)}{\partial \xi}\vert _{\xi=0}\]当我们输入变换后的数据，我们可以写出第k输出值 \(y_k\) 对 \(\xi\) 的导数
\[\frac{\partial y_k}{\partial \xi}\bigg\vert _{\xi=0} = \sum_{i=1}^n \frac{\partial y_k}{\partial x_i}\frac{\partial x_i}{\partial \xi} \bigg\vert_{\xi =0} = \sum_{i=1}^n \frac{\partial y_k}{\partial x_i}\tau_i \bigg\vert_{\xi =0}\]我们就可以用上式作为正则项，即
\[\tilde{L} = L + \lambda \Omega, \qquad where\ \Omega = \sum_n\sum_k(\frac{\partial y_k}{\partial \xi}\bigg\vert _{\xi=0})^2\]其中的 \(\frac{\partial y_k}{\partial x_i}\) 可以通过反向传播得到，\(\tau_i\) 可以通过上面提到的 Finite differences 微分得到。
通过理论可以证明，这两种方法是高度关联的。
在我们之前的模型中，只能为每个输入值预测一个输出值。如果想让模型对每个输入有不同的输出值范围，可以借助混合密度网络(Mixture Density Network, MDN)。MDN 并不是让网络预测单一输出值，而是预测输出的整个概率分布。我们以高斯混合模型来刻画，输出值被建模为多个高斯随机值的和。所以，对于每个输入 x ，我们可以通过该模型得到对应的概率分布函数 \(p(y\vert x)\) ，即
\[p(y\vert x) = \sum_{i=1}^K \pi_k(x)N(y\vert \mathbf{\mu_k(x)},\sigma^2_k(x))\]我们假设某个输入 x 对应的标签值为 K 维向量，在传统的神经网络中，我们的输出层只有 K 个单元，对应标签的维度。而在混合密度网络中
所以混合密度网络最后的输出包含 \(\{a^{\sigma}_k, a^{\pi}_k, a^{\mu}_{kj}\}\) 。而我们对最后的概率分布函数和网络输出之间的对应关系为
\[\pi_k(x) = \frac{exp(a^{\pi}_k)}{\sum_{i=1}^Kexp(a^{\pi}_i)}\] \[\sigma_k(x) = exp(a^{\sigma}_k)\] \[\mu_{kj}(x) = a^{\mu}_{kj}\]接着我们用 w 显式地表示模型中的学习参数，写出混合密度网络中的损失函数
\[L(w) = -\sum_{i=1}^n\ln \{\sum_{k=1}^K \pi_k(x_i,w)N(y\vert \mathbf{\mu_k}(x_i,w),\sigma^2_k(x_i,w))\}\]我们同样可以用误差反向传播来更新参数，这里直接给出理论推导后得到输出层更新公式
\[\frac{\partial L_n}{\partial a^{\pi}_k} = \pi_k - \gamma_k\] \[\frac{\partial L_n}{\partial a^{\mu}_{kj}} = \gamma_k (\frac{\mu_{kj}-y^*_{j}}{\sigma^2_k})\] \[\frac{\partial L_n}{\partial a^{\sigma}_k} = -\gamma_k (\frac{||y^*_n-\mu_k||^2}{\sigma^3_k} - \frac{D}{\sigma_k})\]该方法和 EM 算法的不同是，EM 算法是通过迭代近似得到高斯混合模型中的一系列参数，而混合密度网络是通过神经网络来学习这些参数。