手写数字识别实践（三）：基于BP神经网络与PyQt5

  老师最后给了个用神经网络或者深度学习模型识别手写数字的要求。那就BP神经网络咯（基于误差反向传播算法的前馈神经网络），复习备考期间，一切从简。
  编程及运行环境： VS Code + Win10（Ubuntu也是可以跑的）
  使用软件包： PIL + PyQt5 + numpy + pandas + matplotlib

一、理论基础

1.神经元

  神经元是神经网络的基本构成单元，其典型结构如下所示：

  取 $d$ 维输入样本 $x=(x_1,x_2,...,x_d{)}^T$，就单个神经元而言，其输入和输出的关系为
$$a=f(z)=f(\sum _{i=1}^d{w}_i{x}_i+b)=f(w^{T}x+b)$$   上式中，$a$ 为激活值，$z$ 为净输入，$w=(w_1,w_2,...,w_d{)}^T$ 为权重向量，$b$ 为偏置。

PS：从关系式可以看出，神经元的输入和输出之间并非简单的加权求和，而是在此基础上存在一层映射关系，即激活函数。为什么会有这个东西呢？因为上面展示出来的关系式中，撇开激活函数来看，输入和输出就是线性关系，实际应用中不可能全都是这种理想的线性关系，因此为了提高模型的表达能力，需要引入非线性元素，也就是激活函数。

2.激活函数

  激活函数作为神经元中极其重要的一部分，当然不可能是随便选的。对于一个神经网络而言，所谓的学习其实就是不停地调整各个神经元的输入权重和偏置，使神经网络的输出逐步逼近我们的期望值（PID参数整定既视感有没有！！）。这就要求神经网络具有一种性质：改变某一个权重或偏置很小的值，整个神经网络的输出也应该改变很小（用数学语言来说就是函数连续且可导），否则这种学习就会非常困难。除此之外，函数及其导函数越简单，网络计算效率就会更高；函数的导函数值域也要合适，不能太大也不能太小，否则会影响训练的效率和稳定性。下面学习两类 老师重点强调要复习的 典型的激活函数。

（1）Sigmoid 型激活函数

  Sigmoid 型函数指的是一类 S 型曲线函数（两端饱和函数，左饱和+右饱和）。常用的 Sigmoid 型函数有 Logistic 函数和 Tanh 函数。不过目前网络上的 Sigmoid 函数好像泛指 Logistic 函数。
   Logistic 函数定义为
$$\sigma (x)=\frac{1}{1+e^{-x}}$$   很明显 Logistic 函数是一个单调增函数，其值域为 (0,1)。输入值在0附近时，它可以近似看成线性函数。
  Tanh 函数定义为
$$tanh(x)=\frac{e^x-e^{-x}}{e^x+e^{-x}}=2\sigma (2x)-1$$   从公式就可以看出，Tanh 函数其实就是 Logistic 函数沿 $x$ 轴缩小一倍，沿 $y$ 轴伸展一倍，然后向下平移一格得到的，其值域为 (-1,1)。

PS：Logistic 函数的输出恒大于0。这种输出非零中心化的性质容易使得后一层的神经元输入发生偏置偏移，并进一步减缓梯度下降的收敛速度。另外，如果目的只在于完成手写数字程序设计，只了解 Logistic 函数相关的部分即可。

（2）rectifier 函数（Rectified Linear Unit，ReLU）

  这货实际上是个斜坡函数，定义为
$$ReLU(x)=\{\begin{array}{ll}x,& \text{x >= 0}\\ 0,& \text{x < 0}\end{array}=max(0,x)$$   ReLU 函数为左饱和函数，且在 $x>0$ 时导数恒为1，这在一定程度上环节了神经网络的梯度消失问题，加速梯度下降的收敛速度。但是 ReLU 函数的输出也是非零中心化的，因此和 Logistic 函数存在相同的问题。另外由于死亡 ReLU 问题（即隐藏层的某个 ReLU 神经元可能因为某次不恰当的参数更新而出现在所有训练数据上都不能被激活的情况，这也意味着这个神经元在以后的训练过程中永远都不会被激活），衍生出了下面几个广泛使用的变种 ReLU 。
  Leaky ReLU（带泄露的ReLU）在输入 $x<0$ 时保持一个很小的梯度 $\lambda$ ，这样当神经元在未激活状态也能有一个非零梯度可以更新参数，避免了永久性的“死亡”。函数形式定义为
$$LeakyReLU(x)=\{\begin{array}{ll}x,& \text{x > 0}\\ \gamma x,& \text{x <= 0}\end{array}=max(0,x)+\gamma min(0,x)$$  PReLU（带参数的ReLU）引入一个可学习的参数，不同神经元可以有不同的参数。第 $i$ 个神经元的函数形式定义为
$$PReL{U}_i(x)=\{\begin{array}{ll}x,& \text{x > 0}\\ {\gamma }_{i}x,& \text{x <= 0}\end{array}=max(0,x)+{\gamma }_{i}min(0,x)$$  ELU（指数线性单元）是一个近似零中心化的非线性函数，其定义为
$$ELU(x)=\{\begin{array}{ll}x,& \text{x > 0}\\ \gamma (e^x-1),& \text{x <= 0}\end{array}=max(0,x)+min(0,\gamma (e^x-1))$$  Softplus 函数可以看成 rectifier 函数的平滑版本，其定义为
$$Softplus(x)=log(1+e^x)$$   Softplus 函数的导数刚好是 Logistic 函数。

PS：除去上面的两类激活函数以外，还有 Swish 函数以及 Maxout 单元。

（3）常见激活函数及其导数（公式有用）：

3.网络结构(前馈网络)

  说完单个的神经元，接下来则是多个神经元组成的神经网络。典型的神经网络结构有前馈神经网络、反馈神经网络和竞争学习型神经网络这三种，不过在这里只介绍要用到的BP神经网络，即前馈神经网络结构。在前馈神经网络中，各神经元分属不同的层。每一层的神经元可以接收前一层神经元的信号，并产生信号输出到下一层。第0层叫输入层，最后一层叫输出层，其它中间层叫隐藏层。整个网络中无反馈，信号从输入层到输出层单向传播。下图所示的网络结构就是最简单的三层（隐层只有一层）前馈网络结构。

  前馈神经网络的信息传播过程可由下面的公式表示：
$$a^{(l)}=f_l(z^{(l)})=f_l(W^{(l)}{a}^{(l-1)}+b^{(l)})$$   用文字表述就是把一个特征向量 $x$ 的各个分量按不同权重加权，再加一个常数项偏置，经激活函数生成输入层激活值（隐层节点输入值）。隐层节点的值经过一个激活函数激活后，获得隐层激活值，隐层激活值仿照输入层到隐层的过程，加权再加偏置，获得输出层值，通过一个激活函数得到最后的输出。这就是“前向过程”。
  除了一般的前向过程外，BP神经网络还有一个“后向过程”：将最终输出与我们给定的期望输出比对，得到偏差信息后根据链式法则将其传递至隐层以更新隐层——输出层的权值和偏置，再从隐层按类似方法传递至输入层来更新输入层——隐层的权值和偏置，循环往复地更新各层间连接的参数，直至最终输出与期望输出相符。这里的后向过程就是误差的反向传播（back propagation，BP），这也是BP神经网络的由来。

4.误差反向传播算法

  这里只贴出随机梯度下降的误差反响传播算法、程序设计的关键公式以及编程角度的权值更新步骤。

  算法4.1中第 $l$ 层的误差项 ${\delta }^{(l)}=f_l^{\prime }(z^{(l)})⨀((W^{(l+1)}{)}^T{\delta }^{(l+1)})$ ，即该层的一个神经元激活函数的导数（梯度）与上一层中所有与该神经元相连的神经元误差项的权重和的点乘。
  反映到程序设计的角度描述BP则可以分成两步（不含前向过程）：
  第一步：后向计算每一层的误差项 ${\delta }^{(l)}$ 。对于输出层来说，误差就是 $f(预测值-真实值)$ ；对于其他层则直接套用上面的公式即可。
  第二步：更新参数。直接照着上面给的公式来就好。可以暂时不考虑正则化，此时将上面的 $\lambda$ 置 0 。

PS：在学习过程中发现一般情况下前馈神经网络中神经元的激活函数都会选择 Sigmoid（Logistic）函数。为什么呢？问得好，我他喵也不知道。

5.BP神经网络的特点

• 输入和输出是并行的模拟量；
• 网络的输入输出关系是各层连接的权因子决定，没有固定的算法（奠定了设计网络结构+调参的基调）；
• 权因子是通过学习信号调节的，这样学习越多，网络越聪明；
• 隐含层越多，网络输出精度越高，且个别权因子的损坏不会对网络输出产生大影响。

二、程序设计实践

1.流程图与数据预处理

  直接沿用博文手写数字识别实践（一）：基于朴素贝叶斯分类器与PyQt5中的设计思路与数据的二值化处理结果。

2.神经网络构建

  源码参考第三条老哥的博客，这里只是做了个简单的封装。

import math
import numpy as np
import pandas as pd

feature_len = 784                                                           #一张图片转为行向量是 28 * 28 = 784 个像素点

class BP_network(object):
    def __init__(self, out_num=10, hid_num=36, w1=None, w2=None, hod_offset=None, out_offset=None, inp_lrate=0.25, hid_lrate=0.25, err_th=0.01):
        self.out_num = out_num                                              # 输出节点数
        self.hid_num = hid_num                                              # 隐层节点数(经验公式:sqrt(输入节点数+输出节点数)+0~9)
        self.w1 = 0.2*np.random.random((feature_len, self.hid_num))- 0.1    # 初始化输入层权矩阵(784*36)
        self.w2 = 0.2*np.random.random((self.hid_num, self.out_num))- 0.1   # 初始化隐层权矩阵(36*10)
        self.hid_offset = np.zeros(self.hid_num)                            # 隐层偏置向量
        self.out_offset = np.zeros(self.out_num)                            # 输出层偏置向量
        self.inp_lrate = inp_lrate                                          # 输入层权值学习率
        self.hid_lrate = hid_lrate                                          # 隐层学权值习率
        self.err_th = err_th                                                # 学习误差门限

  输入节点：这波不使用PCA降维，回归到784个原始特征（如需降维自行插入相关环节，改一改输入节点数和隐层节点数即可）；
  输出节点：取一个 $1\times 10$ 的向量，10位中哪一位的激活值最大则将这一位对应的索引位的值作为预测结果，如1000000000对应0，0100000000对应1；
  隐层节点数：根据网上流传的“经验公式”，即 sqrt(输入节点数 + 输出节点数) + 0~9 中任一常数选定为36；
  学习率：其值大小直接影响参数收敛速度，输入层和隐层均取 0.25（粗调发现0.2~0.3的效果比较好）；
  学习误差门限：提前结束训练需要达到的误差精度，取0.01；
  输入层和隐层权矩阵初始化：由于激活函数使用的是 Sigmoid 型中的 Logistics 函数，因此为使模型具有良好的学习特性，将权矩阵初始化在函数中心部分较小的一个范围内（-0.1，0.1）；
  输入层和隐层偏置向量：从0开始学习。

3.神经网络训练

  这里并没有使用学习误差门限提前结束训练，而是跑满6w张图片。每传进一张图就更新一次权值，可以使结果收敛地更好。

	# 获取激活值（激活函数为Logistics函数）
    def get_act(self,x):
        act_vec = []
        for i in x:
            act_vec.append(1/(1+math.exp(-i)))
        act_vec = np.array(act_vec)
        return act_vec

    def get_err(self,e):
        return 0.5*np.dot(e,e)

    # 训练——可使用err_th与get_err() 配合，提前结束训练过程
    def Train(self, trainset, train_labels):
        for count in range(0, len(trainset)):
            t_label = np.zeros(self.out_num)
            t_label[train_labels[count]] = 1

            #前向过程
            hid_value = np.dot(trainset[count], self.w1) + self.hid_offset       # 隐层值
            hid_act = self.get_act(hid_value)                                    # 隐层激活值
            out_value = np.dot(hid_act, self.w2) + self.out_offset               # 输出层值
            out_act = self.get_act(out_value)                                    # 输出层激活值

            #后向过程
            e = t_label - out_act                                                       # 输出值与真值间的误差
            out_delta = e * out_act * (1-out_act)                                   # 输出层delta计算
            hid_delta = hid_act * (1-hid_act) * np.dot(self.w2, out_delta)          # 隐层delta计算
            for i in range(0, self.out_num):
                self.w2[:,i] += self.hid_lrate * out_delta[i] * hid_act             # 更新隐层到输出层权向量
            for i in range(0, self.hid_num):
                self.w1[:,i] += self.inp_lrate * hid_delta[i] * trainset[count]     # 更新输入层到隐层的权向量

            self.out_offset += self.hid_lrate * out_delta                           # 输出层偏置更新
            self.hid_offset += self.inp_lrate * hid_delta

4.神经网络测试

  这部分代码包含统计测试正确率、绘制ROC曲线及计算AUC。这次博主并没有选择使用相关库来计算 TPR 和 FPR，而是直接将神经网络对应输出位的激活值在这一次总输出中的占比视作该测试样本的概率估计值，统计 TP 点和 FP 点。阈值的取法可能是有些问题的（这里的阈值是均匀取的，但实际上概率分布很不均匀），所以绘制出来的 ROC 曲线存疑（自己都不信手撸的模型能有这种性能Orz），在此先挖个必定鸽掉的坑；AUC计算则是使用近似积分的方式估值（纵轴不取平均是因为取一下发现AUC也等于1了QAQ）。

    def Test_Predict(self, testset, test_labels):
        pre_right = np.zeros(10)                    # 每一类预测正确的数量
        act_numbers = np.zeros(10)                  # 每一类的样本实际数量（P）                
        predict = []                                # 预测结果
        accuracy = []                               # 总体预测正确率
        tpr = np.zeros((10,11))              		# 各类在0~1间不同阈值下的tp数，阈值步长0.1
        fpr = np.zeros((10,11))              		# 各类在0~1间不同阈值下的fp数，阈值步长0.1
        auc = np.zeros(10)

        # 统计测试数据中各类数字的数量
        for i in test_labels:
            act_numbers[i] += 1
        n_numbers = np.full(10,10000) - act_numbers      # 不为某一类的样本实际数量（N）

        for count in range(len(testset)):
            hid_value = np.dot(testset[count], self.w1) + self.hid_offset    # 隐层值
            hid_act = self.get_act(hid_value)                                # 隐层激活值
            out_value = np.dot(hid_act, self.w2) + self.out_offset           # 输出层值
            out_act = self.get_act(out_value)                                # 输出层激活值
            predict.append(np.argmax(out_act))

            # 统计正确率
            if np.argmax(out_act) == test_labels[count]:
                pre_right[test_labels[count]] += 1
            if (count+1) % 500 == 0:
                accuracy.append(float(pre_right.sum())/(count+1))            # 计算总体预测准确率
            
            # 每预测完一个测试样本，就统计它在0~1这11个阈值下对于每一类而言分别会归属到tp还是fp
            for i in range(0,10,1):
                for threshold in range(0,11,1):                                    
                    if out_act[i]/sum(out_act) > float(threshold)/10:
                        if int(np.argwhere(out_act == out_act[i])) == test_labels[count]:
                            tpr[i][threshold] += 1/act_numbers[i]
                        else:
                            fpr[i][threshold] += 1/n_numbers[i]
                            
        for i in range(0,10,1):
            for j in range(1,11,1):        
                auc[i] += (fpr[i][j-1] - fpr[i][j]) * tpr[i][j]

        print(pre_right)
        print('-------------------------------------------------------------------')
        print(act_numbers)
        print('-------------------------------------------------------------------')
        print(pre_right/act_numbers)

        plt.figure()
        plt.plot(np.arange(1,21),accuracy,marker="o",markersize=8)
        plt.xlabel("x -- 1:500")
        plt.title("Predict Accuracy Rate")
        plt.show()

        plt.figure()
        for i,color in zip(range(0,10,1),['red','orange','gold','green','cyan','blue','purple','peru','brown','black']):
            plt.plot(fpr[i],tpr[i],color=color,label='Class {0}(area = {1:0.2f})'.format(i,auc[i]))
        plt.plot([0, 1], [0, 1], 'k--', lw=2)
        plt.xlim([0.0, 1.0])
        plt.ylim([0.0, 1.05])
        plt.xlabel('False Positive Rate')
        plt.ylabel('True Positive Rate')
        plt.title('ROC Curve')
        plt.legend(loc="lower right")
        plt.show()

        return accuracy, np.array(predict)

三、实际效果

  7w张图跑满也只花了35s左右，如果再加上学习误差门限的限制还能再砍10s，而识别的准确率基本不受影响。
  未使用误差门限提前结束训练：

  使用误差门限提前结束训练：

  Emmm……这几个模型跑下来都是这个走向。现在我开始怀疑这和测试集有关……

  10类数字识别的ROC曲线（?）绘制情况如下图。可以看到，本文使用的BP神经网络对数字2、5、7、8的分类性能相比其他数字更弱一些，而数字0和1则效果显著（然而自己写的0还是有点拉跨QAQ）。这一点倒是和各类准确率统计反映出来的情况相同。

四、参考资料

• 邱锡鹏《神经网络与深度学习》第4章 前馈神经网络
• BP神经网络识别手写数字项目解析及代码
• 【机器学习】BP神经网络实现手写数字识别
• 深度学习中神经网络的几种权重初始化方法
• ROC原理介绍及利用python实现二分类和多分类的ROC曲线

五、复习阶段追加——经典数值优化算法

  其实在学习模式识别相当长的一段时间内咱都没弄明白为什么会扯到数值优化算法上去（感觉就是脑子懵懵的情况下老师突然就讲了），也是到复习阶段学着学着才有了自己的理解。这玩意儿其实是这么用上的：
  众所周知，设计分类器其实就是设计决策函数。大部分情况下决策函数的参数未知待定，然后就会自然而然地转换到求取最优参数，使模型的性能最优的问题。这时候呢，就可以通过构造一个名为损失函数的东西（故名思义就是出现决策错误时带来的损失的一种量化指标），让要求取的最优参数成为损失函数的极小（大）值点，此时就可以对损失函数使用数值优化算法求解参数啦。

1.梯度下降法

  老师反复强调了这玩意儿的重要性，应该是一个必考点没跑了。搞它！
  梯度下降法是一种常用的一阶优化方法，它的主要目的是通过迭代找到目标函数的最小值，或者收敛到最小值。基本思路还是很简单的，在目标函数连续可微且不知道函数图像具体情况的条件下，从函数上的任意位置出发，只要沿着函数一步步往更低的方向走，最终就能走（收敛）到局部极小点。那么现在就有一个问题：怎么确定哪个方向是“更低”的方向？

  根据泰勒一阶展开式可以得到
$$f(x+\Delta x)\simeq f(x)+\Delta x\nabla f(x)$$  其中 $\Delta x$ 就是走的步长，而 $\nabla f(x)$ 则为目标函数 $f(x)$ 在 $x$ 这一点的梯度（导数）。如果是沿着下降的方向走，那么就应该有
$$f(x+\Delta x)<f(x)$$   由此可知应保证 $\Delta x\nabla f(x)<0$。如果取 $\Delta x=-\alpha \nabla f(x)(\alpha >0)$，则只要 $\nabla f(x)\ne 0$（真要等于0了那就已经完事了），就一定有
$$\Delta x\nabla f(x)=-\alpha (\nabla f(x){)}^2<0$$  因此 $x$ 的更新算式也已经很明确了：
$$x^{\prime }=x-\alpha \nabla f(x)$$

PS：如果令 $\Delta x=\alpha \nabla f(x)(\alpha >0)$，那就是梯度上升法，就是去找目标函数的局部极大点了。另外，对于单变量函数 $f(x)$ 而言，$\nabla f(x)$ 就是该函数对 $x$ 的导数；而对于多变量函数 $f(x_1,x_2)$ 而言，$\nabla f(x)$ 就是该函数分别对 $x_1$ 和 $x_2$ 求偏导后组成的向量。

  根据使用的数据量的不同，梯度下降可以分为三种方法：批量梯度下降（Batch Gradient Descent，BGD）、随机梯度下降（Stochastic Gradient Descent， SGD）和小批量梯度下降（Mini-Batch Gradient Descent, MBGD）。
  批量梯度下降每次都使用训练集中的所有样本来更新参数，因此可以得到全局最优解，但是相对地，训练时间就会比较长；随机梯度下降已经在上面的算法4.1中有介绍，每次只从样本集中随机选择一组数据进行梯度下降，这样训练速度有所提升，不过得到的可能只是局部最优解；小批量梯度下降则算是前两者的一种折中。

2.牛顿迭代法

  牛顿迭代法是一种函数逼近法，其基本思想是在极小点附近用二阶泰勒多项式来近似代替目标函数 $f(x)$，从而求出 $f(x)$ 极小点的估计值。
  根据极值定理，函数取得极值的必要条件是导数（对于多元函数是梯度）为0 ，即
$$f^{\prime }(x)=0$$   而在极值点附近的点 $x_k$处将函数 $f(x)$ 按二阶泰勒展开得到
$$f(x)\simeq f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2$$   两边对 $x$ 求导，忽略三阶项，则
$$0=f^{\prime }(x)\simeq f^{\prime }(x_k)+f^{\prime \prime }(x_k)(x-x_k)$$ $$x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$$   以上就是一元目标函数的牛顿迭代法的迭代算式。将其推广到多元目标函数时，迭代算式就变成了
$$x_{k+1}=x_k-\frac{\nabla f(x_k)}{{\nabla }^{2}f(x_k)}=x_k-\frac{g_k}{H_k}=x_k-H_k^{-1}{g}_k$$  易知 $H_k、g_k$ 分别为函数在 $x_k$ 处的二阶导数和梯度。其中， $H_k$ 也叫 $Hessian$ 矩阵。另外，也不难从推导过程中看出，上述推导过程成立的条件是 $H_k$ 可逆；且只有在 $H_k$ 正定（$H_k>0$）时，上述迭代过程才是下降方向的迭代。

PS：因为 $H_k$ 是一个 $n\ast n$ 的矩阵，每轮迭代都需要求解 $H_k、H_k^{-1}$，这就使得计算的复杂度极大地提高。所以人们期望寻找到能够近似替代 $H_k$ 的矩阵 $G_k$ 来避免大量的矩阵求逆计算，这就是拟牛顿法。