前言
在之前的学习中，已经了解了神经网络及其神经网络的反向传播算法的具体算法分析，现在，我们可以用使用该算法实现手写数字的识别，即也就是用python实现一个多分类问题。
神经网络的实现
首先，根据神经网络的相关算法原理，用python实现神经网络的前向传播算法
- 导入相关库
与之前的代码一样，我们需要导入相关库以便实现一些神经网络的运算，导入以下相关库。
import matplotlib.pyplot as plt
import numpy as np
import scipy.io as scio
import scipy.optimize as opt
- 初始化一些参数
这些参数的介绍如下注释所示：
input_layer_size = 400 # 20x20 输入像素矩阵
hidden_layer_size = 25 # 25个隐藏层
num_labels = 10 # 0到9的十个输出
- 数据可视化
参考python手写数字识别的训练集介绍，导入训练集，并随机选择100个训练集实现数据可视化，具体代码如下所示(参考注释)：
def display_data(x):
(m, n) = x.shape
# m = 100
# n = 400
# 设置每个数字的宽度与高度(像素)
example_width = np.round(np.sqrt(n)).astype(int)# example_width=20
example_height = (n / example_width).astype(int) #example_height=20
# 计算显示的行数与列数
display_rows = np.floor(np.sqrt(m)).astype(int) #display_rows=10
display_cols = np.ceil(m / display_rows).astype(int)#display_rows=10
# 单个图像之间的间隔
pad = 1
# 设置并初始化显示像素矩阵的大小211*211 ,1+(10*20+1)
display_array = - np.ones((pad + display_rows * (example_height + pad),
pad + display_rows * (example_height + pad)))
# 将每个训练样本显示在矩阵中
curr_ex = 0
for j in range(display_rows):
for i in range(display_cols):
if curr_ex > m:
break
# 每次每行和每列读取一个20*20像素的数字，矩阵大小加21
# 实际上矩阵形式可以认为 10*10*400(20*20像素)
# 填充要显示的像素矩阵
max_val = np.max(np.abs(x[curr_ex]))
display_array[pad + j * (example_height + pad) + np.arange(example_height),
pad + i * (example_width + pad) + np.arange(example_width)[:, np.newaxis]] = \
x[curr_ex].reshape((example_height, example_width)) / max_val
curr_ex += 1
if curr_ex > m:
break
# Display image
plt.figure()
plt.imshow(display_array, cmap='gray', extent=[-1, 1, -1, 1])
plt.axis('off')
rand_indices = np.random.permutation(range(m))
selected = X[rand_indices[0:100], :] #随机选择100个数据
display_data(selected)
最后，随机选择100个输入样本，调用以上函数，实现数据可视化，如下所示：
-
神经网络模型选择
我们选择的神经网络如下图所示，由三层构成，包括一个输入层，隐藏层和输出层，输入层是一个数字图片的像素矩阵，大小为20X20像素，所以输入层有400个输入单元（不计算偏置单元），训练集中已经给我们提供了参数矩阵，我们选择的神经网络是一个三层神经网络，显然，需要两个参数矩阵，第一个参数矩阵是25X401维的矩阵，而第二个参数矩阵
是26X10维的矩阵。（与输出层10个数字类别相对应）。
-
反向传播算法的损失函数
神经网络正则化的损失函数如下所示：
很明显，，所以
的输出如下所示：
假设通过神经网络输出的值是5，则
,即向量
的第5个元素是1，其他元素是0.
综上所述，参考反向传播算法损失函数的相关介绍,损失函数与梯度的计算代码如下：
def nn_cost_function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lmd):
theta1 = nn_params[:hidden_layer_size * (input_layer_size + 1)].reshape(hidden_layer_size, input_layer_size + 1)
theta2 = nn_params[hidden_layer_size * (input_layer_size + 1):].reshape(num_labels, hidden_layer_size + 1)
m = y.size
cost = 0
theta1_grad = np.zeros(theta1.shape) # 25 x 401
theta2_grad = np.zeros(theta2.shape) # 10 x 26
Y = np.zeros((m, num_labels)) # 5000 x 10
for i in range(m):
Y[i, y[i]-1] = 1
a1 = np.c_[np.ones(m), X] # 5000 x 401
a2 = np.c_[np.ones(m), sigmoid(np.dot(a1, theta1.T))] # 5000 x 26
hypothesis = sigmoid(np.dot(a2, theta2.T)) # 5000 x 10
reg_theta1 = theta1[:, 1:] # 25 x 400
reg_theta2 = theta2[:, 1:] # 10 x 25
cost = np.sum(-Y * np.log(hypothesis) - np.subtract(1, Y) * np.log(np.subtract(1, hypothesis))) / m \
+ (lmd / (2 * m)) * (np.sum(reg_theta1 * reg_theta1) + np.sum(reg_theta2 * reg_theta2))
e3 = hypothesis - Y # 5000 x 10
e2 = np.dot(e3, theta2) * (a2 * np.subtract(1, a2)) # 5000 x 26
e2 = e2[:, 1:] # drop the intercept column # 5000 x 25
delta1 = np.dot(e2.T, a1) # 25 x 401
delta2 = np.dot(e3.T, a2) # 10 x 26
p1 = (lmd / m) * np.c_[np.zeros(hidden_layer_size), reg_theta1]
p2 = (lmd / m) * np.c_[np.zeros(num_labels), reg_theta2]
theta1_grad = p1 + (delta1 / m)
theta2_grad = p2 + (delta2 / m)
grad = np.concatenate([theta1_grad.flatten(), theta2_grad.flatten()])
return cost, grad
初始化一些参数，并且调用上述函数，计算结果如下所示：
data = scio.loadmat('ex4weights.mat')
theta1 = data['Theta1']
theta2 = data['Theta2']
lmd = 1
nn_params = np.concatenate([theta1.flatten(), theta2.flatten()])
cost, grad = nn_cost_function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lmd)
print(cost)
- sigmoid函数求导
激活函数的求导可以用以下公式表示：
根据以上公式，用代码实现方式如下所示：
def sigmoid_gradient(z):
g = np.zeros(z.shape)
g = sigmoid(z)*(1-sigmoid(z))
return g
随机输入一组参数，计算结果如下所示：
- 随机初始化
正如之前的介绍一样，神经网络中参数初始化与逻辑回归中不同，不能初始化为0，应该是随机初始化。如下代码所示：
def rand_initialization(l_in, l_out):
w = np.zeros((l_out, 1 + l_in))
ep_init = 0.12
w = np.random.rand(l_out, 1 + l_in) * (2 * ep_init) - ep_init
return w
- 梯度检测
为了有效减小神经网络反向传播算法中的造成的误差，为此比较数值化的计算方法与梯度下降算法是否接近，数值化的梯度计算方法如下所示：
def compute_numerial_gradient(cost_func, theta):
numgrad = np.zeros(theta.size)
perturb = np.zeros(theta.size)
e = 1e-4
for p in range(theta.size):
perturb[p] = e
loss1, grad1 = cost_func(theta - perturb)
loss2, grad2 = cost_func(theta + perturb)
numgrad[p] = (loss2 - loss1) / (2 * e)
perturb[p] = 0
return numgrad
在经过梯度检测之后，我们要对权重参数重新进行校正，详细代码如下所示：
def debug_initialize_weights(fan_out, fan_in):
w = np.zeros((fan_out, 1 + fan_in))
w = np.sin(np.arange(w.size)).reshape(w.shape) / 10
return w
有了以上代码，梯度检测算法的实现如下所示：
def check_nn_gradients(lmd):
input_layer_size = 3
hidden_layer_size = 5
num_labels = 3
m = 5
# 产生一些随机参数
theta1 = debug_initialize_weights(hidden_layer_size, input_layer_size)
theta2 = debug_initialize_weights(num_labels, hidden_layer_size)
X = debug_initialize_weights(m, input_layer_size - 1)
y = 1 + np.mod(np.arange(1, m + 1), num_labels)
# 参数展开
nn_params = np.concatenate([theta1.flatten(), theta2.flatten()])
def cost_func(p):
return nn_cost_function(p, input_layer_size, hidden_layer_size, num_labels, X, y, lmd)
cost, grad = cost_func(nn_params)
numgrad = compute_numerial_gradient(cost_func, nn_params)
print(np.c_[grad, numgrad])
以上代码中cost_func()函数和grad_func()函数的实现过程如下所示：
def cost_func(p):
return nn_cost_function(p, input_layer_size, hidden_layer_size, num_labels, X, y, lmd)[0]
def grad_func(p):
return nn_cost_function(p, input_layer_size, hidden_layer_size, num_labels, X, y, lmd)[1]
通过以上代码的运行，得到的反向传播得到的梯度参数和数值化计算得到的梯度参数如下所示,通过分析，可以看出来，两者大小基本一致。
- 查看隐藏层
通过调用display_data()函数，传入隐藏层权重参数，隐藏层的实现效果如下所示：
数据预测
通过输入参数和已经通过神经网络训练之后得到的权重参数，就可以实现手写数字的识别了，具体代码如下,f返回值为训练精确度。
def predict(theta1, theta2, x):
m = x.shape[0]
x = np.c_[np.ones(m), x]
h1 = sigmoid(np.dot(x, theta1.T))
h1 = np.c_[np.ones(h1.shape[0]), h1]
h2 = sigmoid(np.dot(h1, theta2.T))
p = np.argmax(h2, axis=1) + 1
return p
最后，运行以上代码，得到的训练精确度如下所示：
