EX4:手写数字识别（神经网络反向传播）

EX4:手写数字识别（神经网络反向传播）

上一个练习实现了神经网络的前馈传播，并使用它来预测已提供权重的手写数字图像。在本练习中，将实施反向传播算法来学习神经网络的参数。

1.数据读取

和前面的练习类似

import numpy as np
from scipy.io import loadmat

datas=loadmat('ex4data1.mat')
X=datas['X']
y=datas['y']
#为X添加偏置列
new_col=np.ones((X.shape[0],1))
X=np.c_[new_col,X]
#读取参数
weights=loadmat('ex4weights.mat')
Theta1=weights['Theta1']
Theta2=weights['Theta2']
#定义网络参数
hidden_layer_num=25
output_layer_num=10

由于之后用到的优化函数需要输入参数为一维向量，编写merge函数实现参数矩阵的变形及合并。

def merge(param1,param2):
    """ 
    将输入的参数矩阵合并、拉成一维向量
    输入：
        param1:参数矩阵1
        param2:参数矩阵2
    输出：
        params:合并、拉直后参数向量
    """
    shape1=param1.shape
    shape2=param2.shape
    param1=param1.reshape(shape1[0]*shape1[1],order='C')
    param2=param2.reshape(shape2[0]*shape2[1],order='C')
    params=np.append(param1,param2)
    return params

将刚刚读取出的Theta1和Theta2合并成网络参数：

nn_params=merge(Theta1,Theta2)#(10285,)

2.前向传播代价函数

2.1 Sigmoid函数

def sigmoid(z):
    return 1/(1+np.exp(-z))

2.2前馈网络参数计算

改写一下练习3中的神经网络前馈计算函数，返回每一层网络的输入和输出结果。

本次练习中仅采用一个神经元个数为25的隐藏层和神经元个数为10的输出层。

def feedForward(nn_params,X):
    """
    返回前馈神经网络各层的输入和输出结果
    输入：
        nn_params:参数矩阵合并后结果
        X:训练参数集X
    输出：
        a1:输入X
        z2:X经参数矩阵Theta1计算后结果
        a2:z2经sigmoid计算后结果
        z3:a2经参数矩阵Theta2计算后结果
        a3:z3经sigmoid计算后结果
    """
    Theta1=nn_params[:X.shape[1]*hidden_layer_num].reshape(hidden_layer_num,X.shape[1])
    Theta2=nn_params[X.shape[1]*hidden_layer_num:].reshape(output_layer_num,hidden_layer_num+1)
    a1=X
    #hidden layer
    z2=X.dot(Theta1.T)#(5000，25)
    a2=sigmoid(z2)#(5000，25)
    #add bias
    ones_col=np.ones((a2.shape[0],1))
    a2=np.c_[ones_col,a2]#(5000,26)
    #output layer
    z3=a2.dot(Theta2.T)#(5000,10)
    a3=sigmoid(z3)#(5000,10)
    return a1,z2,a2,z3,a3

2.3代价函数

现在将实现神经网络的代价函数。
代价函数：
$$\begin{gathered} J(\theta )=\frac{1}{m}\sum _{i=1}^m\sum _{k=1}^K[-y_k^{(i)}log((h_{\theta }(x^{(i)}){)}_k)-(1-y_k^{(i)})log(1-(h_{\theta }(x^{(i)}){)}_k)] \\ +\frac{\lambda }{2m}[\sum _{j=1}^{25}\sum _{k=1}^{400}({\Theta }_{j,k}^{(1)}{)}^2+\sum _{j=1}^{10}\sum _{k=1}^{25}({\Theta }_{j,k}^{(2)}{)}^2] \end{gathered}$$

参数解释：

• m：训练样本的个数
• K：类别数（最后一层的neuron个数）
• 400：特征个数
• 25：hidden layer中neuron个数（不算偏置）
• 10：output layer中neuron个数

编写代价函数前，首先明确一下代价函数计算的过程：
$$神经网络预测结果：output{[}^{\prime }z{3}^{\prime }]={\left[\begin{array}{cccc}{h}_{\theta }(x^{(1)})& h_{\theta }(x^{(2)})& ...& h_{\theta }(x^{(5000)})\end{array}\right]}^T$$
每一行都是一个10维行向量,shape=(5000,10)
$$编码后的标签集合y：encode\_y=\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ...\\ y^{(5000)}\end{array}\right]$$
每一列都是一个10维行向量,shape=(5000,10)

接下来编写函数实现：

def reEncode(y):
    """
    将标签集y重新编码
    输入：
        y:原标签集，是一个一维列向量
    输出：
        encoded_y:重编码后的标签集，是一个仅含有0和1的标签矩阵
    """
    encoded_y=np.zeros((y.shape[0],10))
    for i in range(y.shape[0]):
        encoded_y[i,y[i,:]-1]=1
    return encoded_y

计算代价函数的第一部分：

我们要做的就是将对应样本$h_{theta}(x^{(i)})$与$y^{(i)}$代入代价函数得到一个结果（结果是一个数值），再将所有的样本得到的结果相加。
如果我们利用矩阵计算就是：
$$-encode\_y\times log(output)-(1-encode\_y)log(1-output)$$
这样我们会得到一个（5000，5000）的矩阵结果，而我们要求encode_y和output对应的样本编号是一样的，所以我们只要该矩阵对角线部分的元素，将这些元素求和并除以样本总数m，就是我们代价函数中第一部分$\frac{1}{m}\sum _{i=1}^m\sum _{k=1}^K[-y_k^{(i)}log((h_{\theta }(x^{(i)}){)}_k)-(1-y_k^{(i)})log(1-(h_{\theta }(x^{(i)}){)}_k)]$的计算值。

计算代价函数的第二部分：

观察代价函数表达式，很容易知道第二部分实际上就是对参数集Theta1,Theta2分别平方以后求和（不算偏置项，即每个参数集的第一列）。

分析完以后，我们编写函数实现代价计算：

def nnCostFunction(nn_params,X,y,reg_param):
    """
    计算代价函数
    输入：
        nn_params:参数向量
        X:训练集
        y:标签集
        reg_param:公式中lambda，用于正则代价函数计算
    输出：
        cost:计算所得代价
    """
    #样本数量与特征数量
    m=X.shape[0]  
    feature_num=X.shape[1]
    #y重编码
    encoded_y=reEncode(y)   
    Theta1=nn_params[:hidden_layer_num*feature_num].reshape(hidden_layer_num,feature_num)
    Theta2=nn_params[hidden_layer_num*feature_num:].reshape(output_layer_num,hidden_layer_num+1)
    #神经网络前馈计算
    output=feedForward(nn_params,X)   
    #代价函数计算
    first=np.diagonal((-encoded_y.dot(np.log(output[-1].T))-(1-encoded_y).dot(np.log(1-output[-1].T))))\
        .sum()/m
    second=reg_param*(np.power(Theta1[:,1:],2).sum()+np.power(Theta2[:,1:],2).sum())/(2*m)
    cost=first+second
    
    return cost 

令reg_param分别取值为0和1，借助给定参数计算一下非正则化和正则化情况下的代价：

J=nnCostFunction(nn_params,X,y,0)#0.2876291651613189
J_regularized=nnCostFunction(nn_params,X,y,1)#0.38376985909092365

输出结果和作业说明中的“非正则化代价0.287629，正则化代价0.383770”相近，可知结果正确，算法编写正确。

3.反向传播

3.1Sigmoid梯度

对Sigmoid函数做梯度计算：
$$\begin{gathered} g(z)=\frac{1}{1+e^{-z}} \\ \frac{dg(z)}{dz}=\frac{e^{-z}}{(1+e^{-z}{)}^2}=g(z)(1-g(z)) \end{gathered}$$
所以Sigmoid的梯度可以复用Sigmoid函数进行计算。

def sigmoidGradient(z):
    return np.multiply(sigmoid(z),1-sigmoid(z))

输入0得到结果为0.25.

3.2 随机初始化

训练神经网络时，随机初始化参数集很重要。

随机初始化的一种有效策略是在范围$[-{ϵ}_{init},{ϵ}_{init}]$均匀地随机选择值。本练习中应当在范围[-0.12,0.12]范围内初始化参数，即选择${ϵ}_{init}=0.12$此值范围可确保参数保持较小并使学习更有效。

选择初始化范围的一个有效策略是将它建立在神经网络的单元数量上，可以令${ϵ}_{init}=\frac{\sqrt{6}}{\sqrt{L_{in}+L_{out}}}$，$L_{in},L_{out}$分别是待初始化参数集的两端神经元数量。

要完成上述的初始化可以使用python方法：np.random.uniform，方法解释如下：

np.random.uniform(low=0.0, high=1.0, size=None)：生成[low,high)内的size个随机数，默认[0,1)之间，size是int型或元组

def randInitializeWeights(L_in,L_out):
	"""
	目标:随机初始化一个便于网络学习的参数集
	输入:
		L_in:参数集左端输入个数
		L_out:参数集右端输出个数
	输出:
		theta:规定范围内随机生成的参数集
	"""
    #范围选择
    epsilon_init=np.sqrt(6)/np.sqrt(L_in+L_out)
    #随机生成矩阵
    theta=np.random.uniform(-epsilon_init,epsilon_init,size=(L_out,L_in+1))
    return theta

随机生成两个参数矩阵，并将它们合并成参数向量：

initial_Theta1=randInitializeWeights(400,25)
initial_Theta2=randInitializeWeights(25,10)
initial_nn_params=merge(initial_Theta1,initial_Theta2)#(10285,)

3.3反向传播算法实现

首先是计算代价函数对于参数的梯度，算法流程如下：

1. 计算前向传播各种参数（a1，z2，a2，z3，a3（即输出h））。
2. 对于output layer的每一个输出单元，设置输出层误差${\delta }_k^{(3)}=(h_k-y_k),y_k\in 0,1$
3. 对于hiddent layer，设置误差为：${\delta }^{(2)}=({\theta }^{(2)}{)}^T{\delta }^{(3)}.\ast g^{,}(z^{(2)})$
4. 计算梯度：$grad=\sum {\delta }^{(l+1)}(a^{(l)}{)}^T$

def backpropagation(nn_params,X,y,reg_param):
    """
    输入：
    	nn_params:网络参数向量
        X:训练数据集(5000,401)
        y:标签集(5000,1)
        reg_param:惩罚参数
    输出：
        J:代价
        grad:梯度展开成的一维向量(长度为401*25+10*26=10285)
    """
    m=X.shape[0]
    feature_num=X.shape[1]
    #y重编码
    encoded_y=reEncode(y)   
    Theta1=nn_params[:hidden_layer_num*feature_num].reshape(hidden_layer_num,feature_num)
    Theta2=nn_params[hidden_layer_num*feature_num:].reshape(output_layer_num,hidden_layer_num+1)
    #前向传播参数计算
    X,z2,a2,z3,a3=feedForward(nn_params,X)#(5000,401),(5000,25), (5000,26), (5000,10), (5000,10)
        
    #计算output layer误差
    delta3=a3-encoded_y#(5000,10)

    #计算hidden layer误差
    delta2=delta3.dot(Theta2[:,1:])*sigmoidGradient(z2)#(5000,25)

    #计算梯度
    grad1=(delta2.T.dot(X)+Theta1)/m#(25,401)
    grad2=(delta3.T.dot(a2)+Theta2)/m#(10,26)
    grad=merge(grad1,grad2)

    #计算代价
    J=nnCostFunction(nn_params,X,y,reg_param)

    return J,grad

4.网络训练

我们借助scipy.optimize中的优化方法minimize进行训练：

from scipy.optimize import minimize
res=minimize(fun=backpropagation,x0=initial_nn_params,args=(X,y,1), 
                method='TNC', jac=True,options={'maxiter':250})

注意到我们设置了迭代次数为250.
根据minimize方法返回的res结果可以得到训练以后的参数向量，对它进行拆分并通过前馈网络计算出预测结果：

params=feedForward(res.x,X)
y_pred=params[-1]#(5000,10)

其中y_pred中每一个元素表示当前行样本识别数值为当前列序号+1的概率，所以我们可以找出每一行中最大的一个数字，它所在的列位置+1就是该行样本的预测结果，可以将该预测结果和标签集进行对比得到训练后神经网络预测的正确率：

y_pred=np.array(np.argmax(y_pred,axis=1)+1)
correct=np.array([1 if y_pred[i]==y[i] else 0 for i in range(y_pred.shape[0])])
accuracy=correct.sum()/y_pred.shape[0]*100
print('Acuuracy:{}%'.format(accuracy))#Acuuracy:99.1%