练习四 - 神经网络

Detection
215 阅读11分钟

练习四 - 神经网络

一、神经网络图示

download.png

二、引言

在本次练习中，我们将会实现一个反向传播算法并将其应用于手写字符识别。同理，我们也需要导入对应的数据文件：ex4data1.mat 与 ex4weights.mat

我们将通过反向传播算法实现神经网络损失函数和梯度计算的非正则化和正则化版本。 我们还将实现随机权重初始化和使用网络进行预测的方法。

三、实现

由于我们在练习3中使用的数据集是相同的，所以我们将重新使用代码来加载数据。

1. 导入数据

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

data = loadmat('ex4data1.mat')
data

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
 '__version__': '1.0',
 '__globals__': [],
 'X': array([[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]),
 'y': array([[10],
        [10],
        [10],
        ...,
        [ 9],
        [ 9],
        [ 9]], dtype=uint8)}

# 查看数据集的维度
data['X'].shape,data['y'].shape

((5000, 400), (5000, 1))

# 提取数据
X = data['X']
y = data['y']

2. 预处理

我们也需要对我们的y标签进行一次one-hot 编码。 one-hot 编码将类标签n（k类）转换为长度为k的向量，其中索引n为“hot”（1），而其余为0。 Scikitlearn有一个内置的实用程序，我们可以使用这个。

如下所示： class 1=[1000000000],class 2=[0100000000],…,class 10=[0000000001]class \ 1 = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix},class \ 2 = \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix},\dots,class \ 10 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{bmatrix}

分类10对应的是手写数字0

from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(sparse=False)
y_onehot = encoder.fit_transform(y)
# 五千个分类向量
y_onehot.shape

(5000, 10)

# 查看One-hot编码前后的第一个标签的维度
y[0], y_onehot[0,:]

(array([10], dtype=uint8), array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]))

# 查看第二个样本的灰度图以及标签
x_0 = X[1].reshape(20, 20)
y_0 = y_onehot[1,:]
plt.imshow(x_0, cmap = 'gray')
plt.xticks([]),plt.yticks([])
plt.show()
y_0

download.png

array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1.])

3. 模型展示

我们的神经网络图示如下：

download.png

一共有三层，第一层是输入层，第二层是隐藏层，第三层是输出层。在本次练习中，我们的输入是一副手写数字灰度图的灰度值阵列。当然，我们导入的X数据已经将20 * 20 灰度图的各灰度值“拉直”为400个输入单元（不包括额外的x_0 = 1）。

我们要为此练习构建的神经网络具有与我们的输入数据（400个像素灰度 + 一个偏置单元 = 401）大小匹配的输入层，25个单位的隐藏层（带有偏置单元的26个），以及一个输出层， 10个单位对应我们的一个one-hot编码类标签向量。

4. Sigmoid函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为：

g(z)=11+e−zg(z) = \frac{1}{1 + e^{-z}}

令z=θTxz = \theta^Tx,则我们的假设函数为：hθ(x)=g(z)=11+e−θTXh_\theta(x) = g(z) = \frac{1}{1+e^{-\theta^TX}}

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

5. 前向传播函数

神经网络的三层单元分布为：

  • 输入层：400 + 1（20 * 20个像素灰度 + 1个偏置单元）
  • 隐藏层：25 + 1
  • 输出层：10（对应热编码的y向量标签）

整个神经网络的图示结构，由此得出前向传播函数：

download.png

def forward_propagate(X, theta1, theta2):
    # 获取样本个数
    m = X.shape[0]
    # 在X中插入全一列
    a1 = np.insert(X, 0, values=np.ones(m), axis=1)
    # 得到隐藏层的输入数据
    z2 = a1 * theta1.T
    # 得到隐藏层的输出数据
    a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)
    # 得到输出层的输入数据
    z3 = a2 * theta2.T
    # 得到输出层的输出数据
    h = sigmoid(z3)
    
    return a1, z2, a2, z3, h

6. 代价函数

回忆可知，不带正则化项的神经网络的代价函数如下：J(θ)=1m∑i=1m∑k=1K[−yk(i)log⁡((hθ(x(i)))k)−(1−yk(i))log⁡(1−(hθ(x(i)))k)]J(\theta)=\frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{K}\left[-y_{k}^{(i)} \log \left(\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)-\left(1-y_{k}^{(i)}\right) \log \left(1-\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)\right]

其中hθ(x(i))h_\theta(x^{(i)})的计算方法就是神经网络的三层，结果就是神经网路输出层的单元个数：10个，并且K = 10，代表了总分类数，而(hθ(x(i)))k=ak(3)(h_\theta(x^{(i)}))_k = a_k^{(3)}代表了对应的向量分量。

我们应该实现实现一个能计算每个样本x(i)x^{(i)}的假设函数值hθ(x(i))h_\theta(x^{(i)})的前馈网络，并求出所有样本的损失函数值之和。

def cost(params, input_size, hidden_size, num_labels, X, y):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 将参数阵列重塑为和神经网络各层维度，因为参数阵列被初始化为行向量
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # 得到前馈神经网络的结果
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # 利用元素乘法简化计算
    J = 0
    for i in range(m):
        # 得到的都是矩阵
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        # 得到标量
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    return J

这个Sigmoid函数我们以前使用过。 前向传播函数计算给定当前参数的每个训练实例的假设函数值。 它的输出形状应该与y的一个one-hot编码相同。

# 初始化设置
input_size = 400
hidden_size = 25
num_labels = 10
learning_rate = 1

# 随机初始化完整网络参数大小的参数数组
params = (np.random.random(size=hidden_size * (input_size + 1) + num_labels * (hidden_size + 1)) - 0.5) * 0.25

# 与代价函数操作一致
m = X.shape[0]
X = np.matrix(X)
y = np.matrix(y)

# 将参数数组解开为每个层的参数矩阵
theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))

theta1.shape, theta2.shape

((25, 401), (10, 26))

# 查看前馈网络的参数情况
a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
a1.shape, z2.shape, a2.shape, z3.shape, h.shape

((5000, 401), (5000, 25), (5000, 26), (5000, 10), (5000, 10))

代价函数在计算假设矩阵h之后，应用代价函数来计算y和h之间的总误差

cost(params, input_size, hidden_size, num_labels, X, y_onehot)

7.469474942991236

7. 正则化代价函数

我们的下一步是增加代价函数的正则化。 它实际上并不像看起来那么复杂 - 事实上，正则化术语只是我们已经计算出的代价的一个补充。 下面是修改后的代价函数：

J(θ)=1m∑i=1m∑k=1K[−yk(i)log⁡((hθ(x(i)))k)−(1−yk(i))log⁡(1−(hθ(x(i)))k)]+λ2m[∑j=125∑k=1400(Θj,k(1))2+∑j=110∑k=125(Θj,k(2))2]\begin{aligned} J(\theta)=& \frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{K}\left[-y_{k}^{(i)} \log \left(\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)-\left(1-y_{k}^{(i)}\right) \log \left(1-\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)\right]+\\ & \frac{\lambda}{2 m}\left[\sum_{j=1}^{25} \sum_{k=1}^{400}\left(\Theta_{j, k}^{(1)}\right)^{2}+\sum_{j=1}^{10} \sum_{k=1}^{25}\left(\Theta_{j, k}^{(2)}\right)^{2}\right] \end{aligned}

def costReg(params, input_size, hidden_size, num_labels, X, y, learningrate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 将参数阵列重塑为和神经网络各层维度，因为参数阵列被初始化为行向量
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # 得到前馈神经网络的结果
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # 利用元素乘法简化计算
    J = 0
    for i in range(m):
        # 得到的都是矩阵
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        # 得到标量
        J += np.sum(first_term - second_term)
    
    J = J / m
    # 计算正则化项值，注意两参数矩阵的起始值，额外添加的一列是不用计算正则化值的
    reg = learningrate /(2 * m) * (np.sum(np.power(theta1[:,1:],2))+np.sum(np.power(theta2[:,1:],2)))
    
    
    return J+reg

costReg(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

7.474770977563316

接下来是反向传播算法。 反向传播会更新参数，这将会减少训练数据上的神经网络误差。 我们需要做的第一件事是计算我们之前创建的Sigmoid函数的梯度的函数。

8. 反向传播

我们的Sigmoid函数的梯度如下：g′(z)=ddzg(z)=g(z)(1−g(z))g^{\prime}(z)=\frac{d}{d z} g(z)=g(z)(1-g(z))

def sigmoid_gradient(z):
    # 由于我们之前使用sigmoid函数的计算结果通常为向量，所以得到的梯度结果也要使用向量乘法
    return np.multiply(sigmoid(z) , (1 - sigmoid(z)))

反向传播神经网络图示：

download.png

  • 反向传播的算法概念： 给定一个样本及标签(x(t),y(t))(x^{(t)},y^{(t)})，我们首先执行前向传播得到整个网络的激活函数，包括了对应样本的假设函数值hΘ(x)h_\Theta(x)，然后，对于第 l 层的第 j 个节点，我们想要计算误差项δj(l)\delta_j^{(l)}，此值用于衡量这个节点对于我们的输出误差大小的影响程度。
  1. 对于输出层，我们直接计算神经网络输出层的激活函数值aj(3)a_j^{(3)} 与对应标签yjy_j的差值作为误差项。
  2. 对于隐藏层，我们会基于第 l + 1层的节点的误差，并将其乘以一个权重作为误差项δj(l)\delta_j^{(l)}

现在我们准备好实现反向传播来计算梯度。 由于反向传播所需的计算是代价函数中所需的计算过程，我们实际上将扩展代价函数以执行反向传播并返回代价和梯度。

  • 反向传播算法的五个步骤：

  1. 对每一个样本x(t)x^{(t)}计算对应的输入层值a(1)a^{(1)}（插入全一列），并执行前向传播，得到第二层、第三层的(z(2),a(2),z(3),a(3))(z^{(2)}, a^{(2)},z^{(3)}, a^{(3)})，注意我们在计算上述值的时候，一定要添加一个全一列来保证a(1),a(2)a^{(1)},a^{(2)}包括误差项。

  2. 对于输出层的每一个单元 k(对应的向量分量) ，我们使用下列式子来计算误差项δk(3)=(ak(3)−yk)\delta_k^{(3)} = (a_k^{(3)}- y_k)

  3. 对于隐藏层(l=2)(l = 2)，我们使用δ(2)=(Θ(2))Tδ(3)⋅∗g′(z(2))\delta^{(2)}=\left(\Theta^{(2)}\right)^{T} \delta^{(3)} \cdot * g^{\prime}\left(z^{(2)}\right)

  4. 使用下列式子计算梯度之和，注意，我们应该跳过δ0(2)\delta_0^{(2)},Δ(l)=Δ(l)+δ(l+1)(a(l))T\Delta^{(l)}=\Delta^{(l)}+\delta^{(l+1)}\left(a^{(l)}\right)^{T}

  5. 通过将累计梯度除以1m\frac{1}{m}得到未正则化的神经网络代价函数的梯度∂∂Θij(l)J(Θ)=Dij(l)=1mΔij(l)\frac{\partial}{\partial \Theta_{i j}^{(l)}} J(\Theta)=D_{i j}^{(l)}=\frac{1}{m} \Delta_{i j}^{(l)}

def backprop(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # initializations
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    
    # perform backpropagation
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        # 计算输出层误差项
        d3t = ht - yt  # (1, 10)
        # 插入全一列，计算隐藏层误差项
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        # 计算累计梯度
        delta1 = delta1 + (d2t[:,1:]).T * a1t
        delta2 = delta2 + d3t.T * a2t
    # 得到最终结果——梯度矩阵    
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    # 将梯度矩阵分解为单个数组
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
    
    return J, grad

反向传播计算的最难的部分（除了理解为什么我们正在做所有这些计算）是获得正确矩阵维度。 顺便说一下，你容易混淆了A * B与np.multiply（A，B）使用。 基本上前者是矩阵乘法，后者是元素乘法（除非A或B是标量值，在这种情况下没关系）。 无论如何，让我们测试一下，以确保函数返回我们期望的。

J, grad = backprop(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)
J, grad.shape

(7.474770977563316, (10285,))

我们还需要对反向传播函数进行一个修改，即将梯度计算加正则化。 最后的正式版本如下。

def backpropReg(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # initializations
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    
    # perform backpropagation
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        # 计算输出层误差项
        d3t = ht - yt  # (1, 10)
        # 插入全一列，计算隐藏层误差项
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        # 计算累计梯度
        delta1 = delta1 + (d2t[:,1:]).T * a1t
        delta2 = delta2 + d3t.T * a2t
    # 得到最终结果——梯度矩阵    
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    # add the gradient regularization term
    delta1[:,1:] = delta1[:,1:] + (theta1[:,1:] * learning_rate) / m
    delta2[:,1:] = delta2[:,1:] + (theta2[:,1:] * learning_rate) / m
    
    # 将梯度矩阵分解为单个数组
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
    
    return J, grad

J, grad = backpropReg(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)
J, grad.shape

(7.474770977563316, (10285,))

我们终于准备好训练我们的网络，并用它进行预测。 这与以往的具有多类逻辑回归的练习大致相似。

9. 应用

from scipy.optimize import minimize

# minimize the objective function
fmin = minimize(fun=backpropReg, x0=params, args=(input_size, hidden_size, num_labels, X, y_onehot, learning_rate), 
                method='TNC', jac=True, options={'maxiter': 250})
fmin

C:\Users\Lonicera_zmy\AppData\Local\Temp\ipykernel_7192\2226448044.py:21: RuntimeWarning: divide by zero encountered in log
  second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))





     fun: 0.33440612969003997
     jac: array([-1.46970727e-04, -2.06746985e-07,  4.29085586e-07, ...,
       -3.45471857e-05, -1.31922585e-04, -2.39683314e-04])
 message: 'Max. number of function evaluations reached'
    nfev: 250
     nit: 20
  status: 3
 success: False
       x: array([ 2.08219381e-01, -1.03373492e-03,  2.14542793e-03, ...,
        3.01993026e-01,  1.20463769e+00,  1.57206779e-01])

由于目标函数不太可能完全收敛，我们对迭代次数进行了限制。 我们的总代价已经下降到0.5以下，这是算法正常工作的一个很好的指标。 让我们使用它发现的参数，并通过网络转发，以获得一些预测。

让我们使用它找到的参数，并通过网络前向传播以获得预测。

X = np.matrix(X)
# 分离一二层的参数矩阵
theta1 = np.matrix(np.reshape(fmin.x[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
theta2 = np.matrix(np.reshape(fmin.x[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
# 得到前向传播的参数
a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
# 得到预测结果
y_pred = np.array(np.argmax(h, axis=1) + 1)
y_pred

array([[10],
       [10],
       [10],
       ...,
       [ 9],
       [ 9],
       [ 9]], dtype=int64)

# 查看维度
y_pred.shape

(5000, 1)

最后，我们可以计算准确度，看看我们训练完毕的神经网络效果怎么样。

# 计算预测精度
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, y)]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))

accuracy = 99.42%

我们已经成功地实现了一个基本的反向传播神经网络，并用它来分类手写数字图像。 在下一个练习中，我们将介绍另一个强大的监督学习算法，支持向量机。

avatar
文章
阅读
粉丝