手写数字识别v1.0

本系列目地：

从0开始写一个神经网路的手写数字识别，随着本系列升级的过程中，加强对基本知识的理解

1. 数据集

方便起见，采用tensorflow的内置数据集，当然也可以在此处下载数据集，

import tensorflow as tf
import matplotlib.pyplot as plt
(x_train, y_train), (x_test, y_test) = tf.keras.datasets.mnist.load_data()
plt.imshow(x_train[0])
print('训练集长度为{}，每一个的大小为{}像素'.format(len(x_train), x_train[0].shape))

结果

训练集长度为60000，每一个的大小为(28, 28)像素

2. 神经网络前向传播

简单的神经网络：

本网络一共三层，第一层可以看做是一个样本的三个属性，记为$a_1^1，a_2^1，a_3^1$，上标代表第一层，下表代表第几个数字，$a$代表本神经元的结果，就是激活值。$a_1^2=\sigma (a_1^1{w}_{11}^2+a_2^1{w}_{12}^2+a_3^1{w}_{11}^2+b_1^2)$，这里$w_{ij}^2$代表着连接第一层到第二层的权重，$b_1^2$代表计算到第二层第一个神经元的偏置，且
$$\sigma (x)=\frac{1}{1+e^{-x}}\text{\hspace{1em}(1)}$$因此，我们可以给出$a_j^{l+1}$:
$$a_j^{l+1}=\sigma (z_j^{l+1})，其中z_j^{l+1}=\sum _k{w}_{kj}^{l+1}{a}_k^l+b_j^{l+1}$$
写成矩阵的形式即为：
$$a^{l+1}=\sigma (z^{l+1}),其中z^{l+1}=w^{l+1}{a}^l+b^{l+1}$$
即，以上图为例子，假设一个样本的维度为$1\ast 3$，通过权重矩阵$3\ast \mathrm{４}$，将其变成一个$1\ast \mathrm{４}$的，然后再通过一个 $5\ast 1$的矩阵，转变为$1\ast 1$

那么假设现在有$n$个样本，则输入为$n\ast 3$，通过第一个权重矩阵计算，就可以有$n$个样本的中间隐藏层，即$n\ast \mathrm{４}$，最后再通过第二个权重矩阵，可以得到$n\ast 1$的一个结果矩阵

在前向传播中，输入是不定的，但是其中的权重矩阵的维度是固定的。其实不难发现，此维度是和神经元的个数是有关系的。

第一个权重矩阵的行的维度肯定等于第一列的神经元个数，列为第二列神经元的个数。
第二个权重矩阵的行的维度肯定等于第二列的神经元个数，列为第三列神经元的个数。
…
第ｎ-1个权重矩阵的行的维度肯定等于第n-1列的神经元个数，列为第n列神经元的个数。

代码：

import numpy as np
class NetWork():
    def __init__(self, layers):
        self.layers = layers
        selfl.weights = [np.random.randn(x, y) for x, y in zip(layers[:-1], layers[1:])]
        self.bias = [np.random.randn(1, y) for y in layers[1:]]
    
    def __call__(self, x):
        for w, b in zip(self.weights, self.bias):
            x = np.dot(x, w) + b
        return x

现在我们给10个样本，建立一个上图的网络

x = np.random.randn(10,3)
model = NetWork([3,5,1])
model(x).shape

可以看到结果即为（10,1）

3. 神经网络反向传播

3.1 负梯度方向

反向传播最主要的目的就是利用梯度下降法来更新这些权重，那么下面讲解为什么可以利用负梯度方向

假设目标函数为$C(v)$，权重参数为$v=(v_1,v_2)$，更新的变化为$\Delta v$, 为了使得$\Delta C=C(v+\Delta v)-C(v)<0$，根据一阶泰勒展开，$$\Delta C\approx C^{\prime }(v_{10})\Delta v_1+C^{\prime }(v_{20})\Delta v_2=\nabla C(v_0)\Delta v$$
这里带一个0的意思表明是一个固定值，不是变量。很明显，当$\Delta v=-\nabla C(v_0)$时，即当更新的变化为负梯度方向的时候，可以有$\Delta C\approx -\nabla C^2<0$，可以保证$C$是减小的。 然后更新$v$的取值
$$v=v+\eta \Delta v=v-\eta \nabla C(v_0)\text{\hspace{1em}(2)}$$这里$\eta$代表着一个学习率。因此在神经网络传播的过程中，我们需要知道梯度的值是多少

3.2 如何求导

求导的过程是反着来的，通过最后一层的导数不停的向前传导。从下面可以看这个过程，最简答的导数，就是对最后一层变量的求导。

3.2.1 标量形式

假设网络有$L$层，则最后一层的输出为
$$a_j^L=\sigma (z_j^L),z_j^L=\sum _k{w}_{kj}^L{a}_k^{L-1}+b_j^L\text{\hspace{1em}(3)}$$
使用最淳朴的目标函数
$$minC=\frac{1}{2n}\sum _j(y_j^{true}-a_j^L{)}^2\text{\hspace{1em}(4)}$$
其向量形式为
$$minC=\frac{1}{2n}||y^{true}-a^L|{|}^2\text{\hspace{1em}(5)}$$
在3.1中，$$\nabla C=(\frac{\partial C}{\partial w_{11}^L},\frac{\partial C}{\partial w_{12}^L},...,\frac{\partial C}{\partial w_{11}^{L-1}},...)\text{\hspace{1em}(6)}$$

为了方便起见，我们记录目标函数对第$L$层第$j$个神经元的导数为
$${\delta }_j^L=\frac{\partial C}{\partial z_j^L}=\sum _k\frac{\partial C}{\partial a_k^L}\frac{\partial a_k^L}{\partial z_j^L}=\frac{\partial C}{\partial a_j^L}{\sigma }^{\prime }(z_j^L)\text{\hspace{1em}(7)}$$
对于最后一层而言，这里$C=C(a_1^L,a_2^L,...,a_n^L)$，$a_j^L=\sigma (z_j^L)$，所以只有当$k=j$的时候，$\partial a_k^L/\partial z_j^L$才不为$0$。根据式子（3），对于最后一层的系数$w_{ij}^L$的导数为：
$$\begin{array}{rl}\frac{\partial C}{\partial w_{ij}^L}& =\frac{\partial C}{\partial z_j^L}\frac{\partial z_j^L}{\partial w_{ij}^L}\\ & ={\delta }_j^L{a}_i^{L-1}\end{array}\text{\hspace{1em}(8)}$$也就是说我们需要保存上一层的输出，来计算对最后一层$w$的求导
同样的根据（3），我们可以得到
$$\begin{array}{rl}\frac{\partial C}{\partial b_j^L}& =\frac{\partial C}{\partial z_j^L}\frac{\partial z_j^L}{\partial b_j^L}\\ & ={\delta }_j^L\end{array}\text{\hspace{1em}(9)}$$
现在我们只要解决对倒数第二层的导数就可以了。因为所有的$L$层神经元都是连接着$L-1$层的任何神经元，因此，对于任何一个倒数第二层的神经元权重，我们有
$$\begin{array}{rlr}\frac{\partial C}{\partial w_{ij}^{L-1}}& =\frac{\partial C}{\partial z_1^L}\frac{\partial z_1^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial w_{ij}^{L-1}}+\frac{\partial C}{\partial z_2^L}\frac{\partial z_2^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial w_{ij}^{L-1}}& +...+\frac{\partial C}{\partial z_k^L}\frac{\partial z_k^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial w_{ij}^{L-1}}\\ & =\sum _k{\delta }_k^L{w}_{kj}^L{\sigma }^{\prime }(z_j^{L-1})a_i^{L-2}& \end{array}\text{\hspace{1em}(10)}$$
这里$$z_j^L=w_{1j}^L{a}_1^{L-1}+w_{2j}^L{a}_2^{L-1}+...+w_{j-1j}^L{a}_{j-1}^{L-1}+w_{jj}^L{a}_j^{L-1}+w_{j+1j}^L{a}_{j+1}^{L-1}+...+w_{kj}^L{a}_k^{L-1}+b_j^L$$$$a_j^{L-1}=\sigma (z_j^{L-1})$$$$z_j^{L-1}=w_{1j}^{L-1}{a}_1^{L-2}+w_{2j}^{L-1}{a}_2^{L-2}+...+w_{i-1j}^{L-1}{a}_{i-1}^{L-2}+w_{ij}^{L-1}{a}_i^{L-2}+w_{i+1j}^{L-1}{a}_{i+1}^{L-2}+...+w_{kj}^{L-1}{a}_k^{L-2}+b_j^{L-1}$$$k$的值是根据上一层的神经元个数确定的。我们还可以把$(10)$重新写为
$$\begin{array}{rl}\frac{\partial C}{\partial w_{ij}^{L-1}}& =\frac{\partial C}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial w_{ij}^{L-1}}\\ & ={\delta }_j^{L-1}{a}_i^{L-2}\end{array}$$因此我们有了每一层的导数递推公式
$${\delta }_j^{l-1}=\sum _k{\delta }_k^l{w}_{kj}^l{\sigma }^{\prime }(z_j^{l-1})\text{\hspace{1em}(11)}$$

以及对$b_j^{L-1}$的导数：
$$\begin{array}{rl}\frac{\partial C}{\partial b_j^{L-1}}& =\frac{\partial C}{\partial z_1^L}\frac{\partial z_1^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial b_j^{L-1}}+\frac{\partial C}{\partial z_2^L}\frac{\partial z_2^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial b_j^{L-1}}+...+\frac{\partial C}{\partial z_k^L}\frac{\partial z_k^L}{\partial a_j^{L-1}}\frac{\partial a_j^{L-1}}{\partial z_j^{L-1}}\frac{\partial z_j^{L-1}}{\partial b_j^{L-1}}\\ & =\sum _k{\delta }_k^L{w}_{kj}^L{\sigma }^{\prime }(z_j^{L-1})\\ & ={\delta }_j^{L-1}\end{array}\text{\hspace{1em}(12)}$$其实通过$(8)-(12)$我们可以发现，其实在计算导数的时候，我们求解的重点就是$\delta$，以及保留各个层是输出$a$。

3.2.2 向量形式

把标量的计算改为向量的计算，可以有效的提高计算的速度。我们接下来把$3.2.1$的标量方程写为向量的形式：
目标函数的向量形式已经由$(5)$给出，对于导数的定义$(7)$式，假设L层有n个神经元，我们可以写为
$$\begin{array}{rl}{\delta }^L& =({\delta }_1^L,{\delta }_2^L,...,{\delta }_{n-1}^L,{\delta }_n^L)\\ & =(\frac{\partial C}{\partial a_1^L}{\sigma }^{\prime }(z_1^L)，\frac{\partial C}{\partial a_2^L}{\sigma }^{\prime }(z_2^L),...,\frac{\partial C}{\partial a_{n-1}^L}{\sigma }^{\prime }(z_{n-1}^L),\frac{\partial C}{\partial a_n^L}{\sigma }^{\prime }(z_n^L))\\ & ={\nabla }_{a^L}C\odot {\sigma }^{\prime }(z^L)\end{array}\text{\hspace{1em}(12)}$$这里$\odot$表示向量相应的位置相乘，${\nabla }_{a^L}$代表对最后一层的梯度。

下面写一下对变量$w,b$的求导的向量书写形式
根据式子$(8)$,　并假设本层的神经元个数为$n$，上一层的神经元个数$m$
$$\begin{array}{rl}\frac{\partial C}{\partial w^L}& =\left(\begin{array}{c}\frac{\partial C}{\partial w_{11}^L}\frac{\partial C}{\partial w_{12}^L}...\frac{\partial C}{\partial w_{1m}^L}\\ \frac{\partial C}{\partial w_{21}^L}\frac{\partial C}{\partial w_{22}^L}...\frac{\partial C}{\partial w_{2m}^L}\\ ...\\ \frac{\partial C}{\partial w_{n1}^L}\frac{\partial C}{\partial w_{n2}^L}...\frac{\partial C}{\partial w_{nm}^L}\end{array}\right)\\ & =(a_1^{L-1},a_2^{L-1},...a_m^{L-1}{)}^T({\delta }_1^L,{\delta }_2^L,...,{\delta }_{n-1}^L,{\delta }_n^L)\\ & =(a^{L-1}{)}^T{\delta }^L\end{array}\text{\hspace{1em}(13)}$$
根据式子$(9)$，我们有
$$\begin{array}{rl}\frac{\partial C}{\partial b^L}& =(\frac{\partial C}{\partial b_1^L},\frac{\partial C}{\partial b_2^L},...,\frac{\partial C}{\partial b_{n-1}^L},\frac{\partial C}{\partial b_n^L})\\ & =({\delta }_1^L,{\delta }_2^L,...,{\delta }_{n-1}^L,{\delta }_n^L)\\ & ={\delta }^L\end{array}\text{\hspace{1em}(14)}$$

对于式子$(11)$的向量形式，我们可以逐项观察，首先对于任意的$l，l=1...L-1,$
$${\delta }^l=({\delta }_1^l,{\delta }_2^l,...,{\delta }_{n-1}^l,{\delta }_n^l)$$
观察第一项：$$\begin{array}{rl}{\delta }_1^l& =({\delta }_1^{l+1}{w}_{11}^{l+1}+{\delta }_2^{l+1}{w}_{21}^{l+1}+...+{\delta }_k^{l+1}{w}_{k1}^{l+1}){\sigma }^{\prime }(z_1^l)\\ & =\left(\begin{array}{c}{\delta }_1^{l+1},{\delta }_2^{l+1},...,{\delta }_k^{l+1}\end{array}\right)\left(\begin{array}{c}{w}_{11}^{l+1}\\ w_{21}^{l+1}\\ ...\\ w_{k1}^{l+1}\end{array}\right){\sigma }^{\prime }(z_1^l)\end{array}$$
观察第二项：$$\begin{array}{rl}{\delta }_2^l& =({\delta }_1^{l+1}{w}_{12}^{l+1}+{\delta }_2^{l+1}{w}_{22}^{l+1}+...+{\delta }_k^{l+1}{w}_{k2}^{l+1}){\sigma }^{\prime }(z_2^l)\\ & =\left(\begin{array}{c}{\delta }_1^{l+1},{\delta }_2^{l+1},...,{\delta }_k^{l+1}\end{array}\right)\left(\begin{array}{c}{w}_{12}^{l+1}\\ w_{22}^{l+1}\\ ...\\ w_{k2}^{l+1}\end{array}\right){\sigma }^{\prime }(z_2^l)\end{array}$$
…

因为之前我们所定义的$w_{ij}$矩阵表示的是前层第i个神经元到本层第j个神经元的权重，而在上式$w$的行数是和$l+1$层的神经元个数一样，因此我们需要进行转置，${\delta }^l$可以简单的写为：
$${\delta }^l={\delta }^{l+1}(w^{l+1}{)}^T\odot {\sigma }^{\prime }(z^l)\text{\hspace{1em}(15)}$$

3.3 神经网络传播过程

1. 输入， 其实就是对于第一层的激活值
2. 前向传播，对每一个$l=2,3,...,L$计算相应的$z^l=a^{l-1}{w}^l+b^l$和$a^l=\sigma (z^l)$
3. 计算各层的导数${\delta }^l$，计算$C_{w^l},C_{b^l}$
4. 更新权重，$w^l=w^l-\eta C_{w^l},b^l=b^l-\eta C_{b^l}$

4. 手写数字最简单代码

结果,经过10 轮反复训练， 最高可以达到86%的准确度

大致过程：
对 每一张图片 都进行$3.3$的过程，更新权重，因此一轮就更新了6万次， 10轮更新了60万次。

class NetWork():
    def __init__(self, layers):
        self.layers = layers
        self.weights = [np.random.randn(x, y) for x, y in zip(layers[:-1], layers[1:])]
        self.bias = [np.random.randn(1, y) for y in layers[1:]]
        self.eta = 0.1
        self.z_record = []
        self.delta = []
        self.activates = []
        self.cws = []
        self.cbs = []
        print('weights[0].shape={},weights[1].shape={}'.format(self.weights[0].shape, self.weights[1].shape))
        print('bias[0].shape={},bias[1].shape={}'.format(self.bias[0].shape, self.bias[1].shape))
    
    # 前向传播
    def __call__(self, x):
        for w, b in zip(self.weights, self.bias):
            x = self.sigmoid(np.dot(x, w) + b)
            
        return x
    
    ## sigmoid 函数定义
    def sigmoid(self, x):
        return 1 / (1 + np.exp(-x))
    
    ## 定义 sigmoid 导函数
    def sigmoid_prime(self, x):
        return self.sigmoid(x) * (1 - self.sigmoid(x))
    
    def train_helper(self, x_train, y_train, epochs, x_test=None, y_test=None):
#         转换为one-hot编码
        y_train = pd.get_dummies(y_train)
        
        
        # 添加epoch
        for epoch in range(epochs):
            for i in range(len(x_train)):
    #         for i in range(200):

                self.train(i, x_train[i], y_train.loc[i].tolist()) 
            print("epoch:{}, acc:{}".format(epoch, self.evaluate(x_test,y_test)))
        
        print("training completed")
    
    def evaluate(self, x_test, y_test):
        pre = [np.argmax(self.__call__(i)) for i in x_test]
        return sum(pre == y_test)
        
        

    
    def train(self,index, x_train, y_train):
        x = x_train[np.newaxis,:]
        z_record = [x] # 添加第一列的值，用于更新第一个w和b
        activates = []
        c_ws = []
        c_bs = []
        
        ## 记录导数
        for w, b in zip(self.weights, self.bias):
            z = np.dot(x, w) + b
            z_record.append(z)
            x = self.sigmoid(z)
            activates.append(x)
        

        # 求出delta导数，最后一层，对应于（12）
        delta_record = [self.cost_derivative(y_train, x) * self.sigmoid_prime(z_record[-1])]
        
        if index % 1000 == 0:
            print('index:{}, MSE loss:{},'.format(index, np.mean((y_train - x) ** 2) ))
            
        # 计算（13）、(14)
        c_w = np.dot(z_record[-2].T, delta_record[0])
        c_b = delta_record[0]
        
        c_ws.append(c_w)
        c_bs.append(c_b)
        
        
        
        # 更新权重
        self.weights[-1] = self.weights[-1] - self.eta * c_w  
        self.bias[-1] = self.bias[-1] - self.eta * c_b

        
        
        
        # 记录其他层的导数delta,计算相应的导数，并更新
        for ly in range(2, len(self.layers)):
            # 计算到神经元的导数,对于于（15）

            this_delta = np.dot(delta_record[0], self.weights[-ly+1].T) * self.sigmoid_prime(z_record[-ly]) 
            delta_record.insert(0, this_delta) #
            
            # 计算对本层ｗ　和　ｂ的导数，并且更新
            c_w = np.dot(z_record[-(ly+1)].T, delta_record[0]) # (13)
            c_b = delta_record[0] # (14)
            
            # 更新权重
            self.weights[-ly] = self.weights[-ly] - self.eta * c_w  
            self.bias[-ly] = self.bias[-ly] - self.eta * c_b
            
            c_ws.append(c_w)
            c_bs.append(c_b)
        
        self.z_record = z_record
        self.delta = delta_record
        self.activates = activates
        self.cws = c_ws
        self.cbs = c_bs
        
    ## 定义对目标函数的导数
    def cost_derivative(self, y_true, y_pre):
        return y_pre - y_true    
        
    

        
    

    

        
    
import tensorflow as tf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
(x_train, y_train), (x_test, y_test) = tf.keras.datasets.mnist.load_data()
## plt.imshow(x_train[0])
x_train = x_train.reshape(x_train.shape[0],-1)
x_test = x_test.reshape(x_test.shape[0],-1)
# x_test = (np.max(x_test) - x_test) / np.max(x_test)
print('训练集长度为{}，每一个的大小为{}像素'.format(len(x_train), x_train[0].shape))

model = NetWork([784,13,10])
# 由于数据太大了，很有可能算出来的值太大了，计算的时候容易出现梯度消失啊
# 因此把数据进行归一化
x_train = (np.max(x_train) - x_train) / np.max(x_train)
# 开始训练
model.train_helper(x_train, y_train, 10, x_test, y_test)

结果
weights[0].shape=(784, 13),weights[1].shape=(13, 10)
bias[0].shape=(1, 13),bias[1].shape=(1, 10)
epoch:0, acc:8078
epoch:1, acc:8266
epoch:2, acc:8611
epoch:3, acc:8034
epoch:4, acc:8413
epoch:5, acc:8069
epoch:6, acc:8573
epoch:7, acc:8682
epoch:8, acc:8258
epoch:9, acc:8283