反向传播代码实现

数据可视化

这里的数据集使用的是用矩阵表示的灰度数字图像，其包含两个部分，灰度图像部分$X$和结果部分$y$。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat

path='E:\Coursera-ML-AndrewNg-Notes\code\ex4-NN back propagation\ex4data1.mat'
data=loadmat(path)
X=data['X']
y=data['y']
y=y.astype(np.float64)
X,y

结果如下图

显示灰度图像

def displayData(x):
    fig,ax=plt.subplots(nrows=10,ncols=10,sharey=True,sharex=True)
    # 返回一个规模为(100,),大小在[0,5000)的一维向量
    pick=np.random.randint(0,5000,(100,))
    for row in range(10):
        for col in range(10):
            # matshow即显示矩阵代表的图像,gray_r即灰度图像
            ax[row,col].matshow(x[pick[10*row+col]][:].reshape(20,20).T,cmap='gray_r')
    plt.xticks([]) # 设置x坐标
    plt.yticks([]) # 设置y坐标
    plt.show()
displayData(X)

结果如下图

获取权重矩阵

path1='E:\Coursera-ML-AndrewNg-Notes\code\ex4-NN back propagation\ex4weights.mat'
thetas=loadmat(path1)
#第一层的权重矩阵
theta1=thetas['Theta1']
#第二层的权重矩阵
theta2=thetas['Theta2']
theta1 #显示第一层的权重矩阵

结果如下图

预处理$X$和$y$

即在$X$中添加一列全为$1$的列，当作$x_{0}$，将每个$y$处理为一个$(1,10)$的向量，全部组合起来即为一个$(5000,10)$的矩阵

def init(y):
    m=len(y)
    res=np.zeros((5000,10))
    for i in range(m):
        res[i,int(y[i])-1]=1
    return res
X=np.c_[np.ones(X.shape[0]),X]

$sigmoid$函数实现

\[f(z)=\frac{1}{1+e^{-z}} \]

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

前向传播函数实现

\[a^{(2)}=f(z^{(2)})=f(X\theta^{(1)}) \]

对$a^{(2)}$添加一列$1$，作为$a^{(2)}_{0}$，得到一个新的$a^{(2)}$

\[z^{(3)}=a^{(2)}\theta^{(2)} \]

\[a^{(3)}=f(z^{(3)}) \]

$a^{(3)}$即为我们的预测值，具体实现时需要注意它们的维度

def forward(X,theta1,theta2):
    z2=X@theta1.T #(5000,25)
    a2=sigmoid(z2) 
    a2=np.c_[np.ones(a2.shape[0]),a2] #(5000,26)
    z3=a2@theta2.T #(5000,10)
    a3=sigmoid(z3) #(5000,10)
    return z2,a2,z3,a3

损失函数的实现

未正则化的损失函数

\[J(\theta)=-\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{K}[y^{(i)}_{k}\ln h_{\theta}(x^{(i)})_{k}+(1-y^{(i)}_{k})\ln(1-h_{\theta}(x^{(i)})_{k})] \]

其中$K$为标签总数(输出层的神经元个数)，在这里为$10$，$h_{\theta}(x)即为a^{(3)}$，$h_{\theta}(x^{(i)})$即为输入第$i$个样本的值得到的$a^{(3)}$值

def cost(X,y,theta1,theta2):
    m=len(y)
    y_1=init(y) #(5000,10)
    z2,a2,z3,a3=forward(X,theta1,theta2)
# 最后的损失函数中的y*log(h(x))不是把y_1矩阵和a2直接相乘，
# 而是y_1的第i行与a3的第i行相乘，如果直接相乘，会导致y_1的第i行和a3的所有行都相乘
    error=0
    for i in range(m):
        first=-np.multiply(y_1[i,:],np.log(a3[i,:]))
        second=np.multiply((1-y_1[i,:]),np.log(1-a3[i,:]))
        error+=np.sum(first-second)
    return error/m

对函数进行测试

cost(X,y,theta1,theta2)

结果如下

0.2876291651613187

正则化的损失函数

\[J(\theta)=-\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{K}[y^{(i)}_{k}\ln h_{\theta}(x^{(i)})_{k}+(1-y^{(i)}_{k})\ln(1-h_{\theta}(x^{(i)})_{k})] \]

\[+\frac{\lambda}{2m}[\sum_{i=1}^{25}\sum_{j=1}^{400}(\theta^{(1)}_{ij})^2+\sum_{i=1}^{10}\sum_{j=1}^{25}(\theta^{(2)}_{ij})^2] \]

式中的$\lambda$即代码中的$learningRate$

def cost_reg(X,y,theta1,theta2,learningRate):
    m=len(y)
    error=cost(X,y,theta1,theta2)
    reg=0
    for i in range(len(theta1)):
        reg+=np.sum(np.multiply(theta1[i,1:],theta1[i,1:]))
    for i in range(len(theta2)):
        reg+=np.sum(np.multiply(theta2[i,1:],theta2[i,1:]))
    reg=learningRate*reg/(2*m)
    return error+reg

对函数进行测试

cost_reg(X,y,theta1,theta2,1)

结果如下

0.3837698590909234

反向传播

$sigmoid$梯度函数的实现

\[f^{'}(z)=f(z)(1-f(z)) \]

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

反向传播函数的实现

\[\frac{\partial J(\theta )}{\partial \theta ^{(k)}_{ij}}=\left\{\begin{matrix} \frac{1}{m}\sum_{\mu=1}^{m}[\Delta^{(k)}_{ij}]^{(\mu)}& if \quad j=0\\ & \\ \frac{1}{m}\sum_{\mu=1}^{m}[\Delta^{(k)}_{ij}]^{(\mu)}+\frac{\lambda}{m}\theta_{ij}^{(k)}& if \quad j\ne0 \end{matrix}\right.\]

式子中$m$为样本个数，$[\Delta^{(k)}_{ij}]^{(\mu)}$为第$\mu$个样本在第$k$个权重矩阵中的第$i$行第$j$列的梯度值，即$[\Delta^{(k)}_{ij}]^{(\mu)}=[\delta^{(k+1)}_{i}\cdot a_{j}^{(k)}]^{(\mu)}$，其中

\[\delta^{(k)}_{i}=\left\{\begin{matrix} a^{(k)}_{i}-y_{i}& if \quad k=l\\ & \\ \delta^{(k+1)}\cdot\theta^{(k)}_{\cdot i}\cdot f^{'}(z^{(k)}_{i}) & if \quad 2\le k<l \end{matrix}\right.\]

式子中的$\theta^{(k)}_{\cdot i}$表示第$k$个权重矩阵中的第$i$列，$i$的取值范围为$[0$，第$k$层神经元的个数$]$
具体实现时，可以使用向量来加快计算，首先对于每个样本，单独计算出它们的$\frac{\partial J(\theta )}{\partial \theta ^{(k)}}$值，而后进行求和并取平均，这样便完成了对于$\frac{1}{m}\sum_{\mu=1}^{m}[\Delta^{(k)}]^{(\mu)}$的计算，而后再对正则项进行处理
代码中的$delta2\_i$和$delta3\_i$即第$i$个样本的$\delta^{(2)}$和$\delta^{(3)}$值，$grad1,grad2$即所有样本的$\Delta^{(1)},\Delta^{(2)}$值的和，$\delta^{(2)}=\delta^{(3)}\theta^{(2)}f^{'}(z^{(2)})$，$\delta^{(3)}=a^{(3)}-y$

def back_propagation(params,X,y,learningRate,hidden_size,input_size,num_labels):
    theta1=params[:hidden_size*(input_size+1)]
    theta2=params[hidden_size*(input_size+1):]
    theta1=np.reshape(theta1,(hidden_size,input_size+1))
    theta2=np.reshape(theta2,(num_labels,hidden_size+1))
    # (5000,25) (5000,26)(5000,10) (5000,10) 
    z2,a2,z3,a3=forward(X,theta1,theta2) 
    m=X.shape[0] #5000
    J=cost_reg(X,y,theta1,theta2,learningRate)
    y=init(y) #(5000,10)
    grad1=np.zeros(theta1.shape) # 第一个权重矩阵的梯度 (25,401)
    grad2=np.zeros(theta2.shape) # 第二个权重矩阵的梯度 (10,26)
    for i in range(m):
        x_i=np.mat(X[i,:]) #(1,401)
        z2_i=np.mat(z2[i,:]) #(1,25)
        a2_i=np.mat(a2[i,:]) #(1,26)
        a3_i=np.mat(a3[i,:]) #(1,10)
        y_i=np.mat(y[i,:]) #(1,10)
        delta3_i=a3_i-y_i #第三层的误差 (1,10)
#         print(z2_i.shape,x_i.shape,a2_i.shape,a3_i.shape,y_i.shape)
        z2_i=np.c_[np.ones(1),z2_i] #(1,26)
        #第二层的误差 (1,26)
        delta2_i=np.multiply((theta2.T@delta3_i.T).T,sigmoid_grad(z2_i)) 
        grad1=grad1+(delta2_i[:,1:]).T*x_i
        grad2=grad2+delta3_i.T*a2_i
    grad1=grad1/m
    grad2=grad2/m #对加和求平均
    grad1[:,1:]=grad1[:,1:]+(theta1[:,1:]*learningRate)/m #正则项的处理
    grad2[:,1:]=grad2[:,1:]+(theta2[:,1:]*learningRate)/m
    grad=np.concatenate((np.ravel(grad1),np.ravel(grad2)))
    return J,grad

进行迭代

\[\theta^{(k)}=\theta^{(k)}-\alpha\frac{\partial J(\theta )}{\partial \theta ^{(k)}} \]

iterations=10 #迭代次数
hidden_size=25
input_size=400
num_labels=10
alpha=2 #学习率
params=np.random.random(size=hidden_size*(input_size+1)+num_labels*(hidden_size+1))
params=(params-0.5)*0.24
learningRate=1 # λ值
for i in range(iterations):
    J,grad=back_propagation(params,X,y,learningRate,hidden_size,input_size,num_labels)
    params-=alpha*grad
    print("迭代次数：",i,"损失函数值：",J)

结果如下图

使用工具库来计算最优参数

from scipy.optimize import minimize
learningRate=1
# fun为需要最小化的函数，x0为需要优化的参数，args为最小化函数里面的其他参数
# 且定义fun时，x0代表的参数必须为fun函数中的第一个参数，且shape必须为(n，)
# jac: 计算梯度的函数，且返回的值必须与x0的形状相同
fmin = minimize(fun=back_propagation,x0=(params),args=(X,y,learningRate,hidden_size,input_size,num_labels),method='TNC',jac=True,options={'maxiter':250})
fmin

结果如下图

图中的$x$即为我们所求的最优参数

验证最优参数

损失函数值

tmp1=np.reshape(fmin.x[:hidden_size*(input_size+1)],(hidden_size,input_size+1))
tmp2=np.reshape(fmin.x[hidden_size*(input_size+1):],(num_labels,hidden_size+1))
cost_reg(X,y,tmp1,tmp2,20)

结果如下

3.9666331849091447

最优参数得到的损失函数值比上面迭代得到的损失函数值还要大的原因，以下面这个例子说明
eg：现在来预测图片的猫和狗的类别，若类别为狗，则值为$1$，否则为$0$。设一个样本输出的是$[0.9，0.1]$，即预测猫的概率为$0.1$，狗的概率为$0.9$，另一个样本输出的是$[0.6，0.4]$，假设这两个样本的真实值都是狗，即$y$向量化后均为$[1,0]$。则对于第一个样本来说，其损失函数值$$J(\theta)=-(1,0)\cdot(\ln0.9,\ln0.1)-(1-1,1-0)\cdot(\ln(1-0.9),\ln(1-0.1))$$

\[=-2\ln0.9=0.2107 \]

而对第二个样本来说，其损失函数值

\[J(\theta)=-2*\ln(0.6)=1.021 \]

可以看出这两个样本都会被预测为狗，但其损失函数差别很大，因此在训练时，是有可能会出现损失函数变大的同时准确率也变大的情况的。比如现在有三个二分类的样本，它们的真实值为$1,1,0$，即$y$向量化后为$[1,0],[1,0],[0,1]$，若样本为狗，则值为$1$，反之为$0$。对于第一组参数，它对三个样本预测的概率分别为$[0.99,0.01],[0.99,0.01],[0.501,0.499]$，即它预测第一个样本为狗的概率为$0.99$，猫的概率为$0.01$，后面的以此类推，而对于第二组参数，它对三个样本预测的概率分别为$[0.501,0.499],[0.501,0.499],[0.499,0.501]$。则对于第一组参数其损失函数值和准确率分别为：

\[J(\theta)=-4\ln(0.99)-2\ln(1-0.501)=1.43 \]

\[accuracy=\frac{2}{3}=66.6\% \]

而对于第二组参数，其损失函数值和准确率分别为：

\[J(\theta)=-6\ln(0.501)=4.14 \]

\[accuracy=100\% \]

可以看到对于第二组参数，其损失函数值和准确率均提高了，因此，损失函数值越小，准确率不一定就越高

准确率

z2,a2,z3,a3=forward(X,tmp1,tmp2)
def predict(a3):
    return np.argmax(a3,axis=1)+1
y_pred=predict(a3)
correct=0
for (a,b) in zip(y_pred,y):
    correct+=1 if a==b else 0
correct/5000

结果如下

0.9904

预测值与实际值比较

from sklearn.metrics import classification_report#评价报告包
print(classification_report(y, y_pred))

结果如下图

隐藏层可视化

fig,ax=plt.subplots(nrows=5,ncols=5,sharex=True,sharey=True,figsize=(6,6))
for r in range(5):
    for c in range(5):
        ax[r,c].matshow(np.array(tmp1[5*r+c][1:].reshape((20,20))).T,cmap='gray_r')
plt.xticks([])
plt.yticks([])
plt.show()

结果如下图

posted @ 2022-08-26 21:43 kris-phl 阅读(42) 评论(0) 收藏举报

刷新页面返回顶部

kris-919

反向传播代码实现

数据可视化

显示灰度图像

获取权重矩阵

预处理\(X\)和\(y\)

\(sigmoid\)函数实现

前向传播函数实现

损失函数的实现

反向传播

使用工具库来计算最优参数

预测值与实际值比较

隐藏层可视化

公告