【机器学习练习 4】 - 神经网络反向传播练习题
本章代码涵盖了基于Python的解决方案，用于Coursera机器学习课程的第四个编程练习。
数据集链接
对于这个练习，我们将再次处理手写数字数据集，这次使用反向传播的前馈神经网络。 我们将通过反向传播算法实现神经网络成本函数和梯度计算的非正则化和正则化版本。 我们还将实现随机权重初始化和使用网络进行预测的方法。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
data = loadmat('ex4data1.mat')
data
原始数据是利用matlab储存的，我们也可以转换成csv格式等
# 1. 加载 .mat 文件
data = loadmat('ex4data1.mat') # 确保文件路径正确
# 2. 提取变量 X 和 y
X = data['X'] # 特征矩阵
y = data['y'] # 标签向量
# 3. 处理标签（MATLAB的10表示0）
y = y.flatten() # 从 (5000,1) 转为 (5000,)
y[y == 10] = 0 # 将标签10改为0
# 4. 转换为 DataFrame
df = pd.DataFrame(X, columns=[f'feature_{i}' for i in range(X.shape[1])])
df['label'] = y # 添加标签列
# 5. 保存为 CSV
df.to_csv('data.csv', index=False)
print("CSV 文件已保存！")
CSV 文件已保存！
由于我们以后需要这些（并将经常使用它们），我们先来创建一些有用的变量。
X = data['X']
y = data['y']
X.shape, y.shape#看下维度
((5000, 400), (5000, 1))
我们也需要对我们的y标签进行一次one-hot 编码。 one-hot 编码将类标签n（k类）转换为长度为k的向量，其中索引n为“hot”（1），而其余为0。 Scikitlearn有一个内置的实用程序，我们可以使用这个。
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(sparse_output=False)
y_onehot = encoder.fit_transform(y.reshape(-1, 1)) # 输入需为2D
print(y_onehot.shape) # 应为 (5000, 10)（假设有10个类别）
(5000, 10)
y[0], y_onehot[0,:] #第10个位置为1，其余为0
标签0的第10个位置为1，其余为0
(array([10], dtype=uint8), array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]))
y[4999], y_onehot[4999,:]
标签4999的第9个位置为1，其余为0
(array([9], dtype=uint8), array([0., 0., 0., 0., 0., 0., 0., 0., 1., 0.]))
我们要为此练习构建的神经网络具有与我们的实例数据（400 +偏置单元）大小匹配的输入层，25个单位的隐藏层（带有偏置单元的26个），以及一个输出层， 10个单位对应我们的一个one-hot编码类标签。 有关网络架构的更多详细信息和图像，请参阅“练习”文件夹中的PDF。
我们需要实现的第一件是评估一组给定的网络参数的损失的代价函数。 源函数在练习文本中（看起来很吓人）。 以下是代价函数的代码。
sigmoid 函数
g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为：
g
(
z
)
=
1
1
+
e
−
z
g\left( z \right)=\frac{1}{1+{{e}^{-z}}}
g(z)=1+e−z1
合起来，我们得到逻辑回归模型的假设函数：
h
θ
(
x
)
=
1
1
+
e
−
θ
T
X
{{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-{{\theta }^{T}}X}}}
hθ(x)=1+e−θTX1
def sigmoid(z):
return 1 / (1 + np.exp(-z))
前向传播函数
(400 + 1) -> (25 + 1) -> (10)
def forward_propagate(X, theta1, theta2):
m = X.shape[0] # 样本数
a1 = np.insert(X, 0, values=np.ones(m), axis=1) # (5000, 401)
z2 = a1 * theta1.T # (5000, 401) * (401, 26) -> (5000, 26)
a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1) # (5000, 26)
z3 = a2 * theta2.T # (5000, 26) * (26, 10) -> (5000, 10)
h = sigmoid(z3) # (5000, 10)
return a1, z2, a2, z3, h
代价函数
def cost(params, input_size, hidden_size, num_labels, X, y, learning_rate):
m = X.shape[0]
X = np.matrix(X) # 5000x400
y = np.matrix(y) # 5000x10
# 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)))) # 25x401
theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1)))) # 10x26
# run the feed-forward pass
a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2) # 5000x401 5000x25 5000x26 5000x10 5000x10
# compute the cost
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 # 25*401+10*26
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)))) # 25x401
theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1)))) # 10x26
theta1.shape, theta2.shape
((25, 401), (10, 26))
a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2) # 5000x401 5000x25 5000x26 5000x10 5000x10
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, learning_rate)
6.90818002768805
正则化代价函数
我们的下一步是增加代价函数的正则化。 它实际上并不像看起来那么复杂 - 事实上，正则化术语只是我们已经计算出的代价的一个补充。 下面是修改后的代价函数。
def cost(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)
# compute the cost
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
# 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)))
return J
cost(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)
6.913437587821537
接下来是反向传播算法。 反向传播参数更新计算将减少训练数据上的网络误差。 我们需要的第一件事是计算我们之前创建的Sigmoid函数的梯度的函数。
def sigmoid_gradient(z):
return np.multiply(sigmoid(z), (1 - sigmoid(z)))
现在我们准备好实施反向传播来计算梯度。 由于反向传播所需的计算是代价函数中所需的计算过程，我们实际上将扩展代价函数以执行反向传播并返回代价和梯度。
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)
#通过 forward_propagate 函数计算各层的激活值（a1, z2, a2, z3, h），其中 h 是输出层的预测值。
# 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,:])) # y=1 时的代价函数
second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:])) # y=0 时的代价函数
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) 输出层误差 d3t：
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
# unravel the gradient matrices into a single array
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
(6.913437587821537, (10285,))
我们还需要对反向传播函数进行一个修改，即将梯度计算加正则化。 最后的正式版本如下。
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)
# d2t = np.multiply((theta2.T @ d3t.T).T, sigmoid_gradient(z2t))
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
# unravel the gradient matrices into a single array
grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
return J, grad
J, grad = backprop(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)
J, grad.shape
(6.913437587821537, (10285,))
我们终于准备好训练我们的网络，并用它进行预测。 这与以往的具有多类逻辑回归的练习大致相似。
from scipy.optimize import minimize
# minimize the objective function
fmin = minimize(
fun=backprop,
x0=params,
args=(input_size, hidden_size, num_labels, X, y_onehot, learning_rate),
method='TNC',
jac=True, # 确保 backprop 返回梯度 (J, grad)
options={'maxfun': 250} # TNC 用 maxfun 而非 maxiter
)
fmin
message: Max. number of function evaluations reached
success: False
status: 3
fun: 0.3295327788794764
x: [-1.720e-01 -3.282e-03 ... 6.581e-01 -2.870e+00]
nit: 21
jac: [ 7.187e-06 -6.565e-07 ... 7.333e-05 5.506e-05]
nfev: 250
由于目标函数不太可能完全收敛，我们对迭代次数进行了限制。 我们的总代价已经下降到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)
最后，我们可以计算准确度，看看我们训练完毕的神经网络效果怎么样。
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.53999999999999%
我们已经成功地实施了一个基本的反向传播神经网络，并用它来分类手写数字图像。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
