本文通过BP算法实现MNIST手写数字识别，详细介绍BP算法原理及Python实现过程，包括权重与偏置初始化、前向传播、随机梯度下降等关键步骤。
MNIST手写识别系列——BP算法
这次想使用的是BP算法来进行求解。
关于BP算法的求解，我强烈推荐这篇博文曾梓华——一文详解神经网络 BP 算法原理及 Python 实现
详细推导可以仔细看博文，我这里只是列出几个重要的公式，下面代码讲解会用到。
这个公式是神经元输出值a的计算公式。其中w是权重，b是偏置，f函数采用的是sigmoid激活函数。
这个是均方误差，我们逆向反馈的目标就是让E函数的值尽可能小。接下来会对该值来求偏导。
权重和偏置的更新公式。
权重w的更新量（需要有一系列的推导）。
偏置b的更新量。
δ 值的更新式子。
sigmoid函数求导。
说完了数学的式子（具体的推导建议去看上面推荐的博客），接下来我们来看具体代码实现。
我使用的是anaconda+jupyter notebook的环境。用Keras框架，Keras可以很容易加载MNIST数据集。这里的代码参考：https://github.com/edvardHua/Articles 但是由于作者使用的是外部数据导入MNIST，并且使用的是python2.7，所以我进行了一些修改。使得更加符合自己以往的学习。
需要导入的包
import pdb
import random
import numpy as np
from keras.datasets import mnist
接下来的内容包含在一个类中。但为了方便说明，我把类的各个方法拆开，方便解释。
先来看看整体：
class Network(object):
def __init__(self, sizes):
# 作用：初始化w和b
# param sizes: list 类型，存储每层神经网络的神经元数目
# 例如sizes=[2,3,2]表示输出层有两个神经元，隐藏层有3个神经元，输出层有2个神经元
# num_layers代表神经网络数目
self.num_layers = len(sizes)
self.sizes = sizes
# 初始输入层，随机产生每层中的y个神经元的biase值（0,1）
self.biases = [np.random.randn(y, 1) for y in sizes[1:]]
self.weights = [np.random.randn(y,x) for x, y in zip(sizes[:-1],sizes[1:])]
# biases和weights是共用的参数。只有一套
# biases和weights其实是一个不规则的二维数组。如果输入是3层，那么数组有2行
# 第一行30个维度，第二行10个维度
# 这个weights的初始方法算是学到了。假设sizes是[a1,a2,a3,a4]
# 那么输出是（a1,a2）(a2,a3) (a3,a4)
def feedforward(self, x):
# 作用：向前传输计算每个神经元的值
# param a：输入值
# return：计算后每个神经元的值
a = x.reshape(x.shape[0],1) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
for b, w in zip(self.biases, self.weights):
# 加权求和再加上biase
a = sigmoid(np.dot(w,a)+b)
# dot是矩阵乘法的意思
return a
def SGD(self, x_train, y_train, x_test, y_test, epochs, mini_batch_size, eta):
# 作用：随机梯度下降
# param training_data：输入的训练集
# param epochs：迭代次数
# param mini_batch_size：小样本数量
# param eta：学习率
# param test_data：测试数据集
n_test = x_test.shape[0]
n_train = x_train.shape[0] # n_train是总的训练样本数目
for j in range(epochs):
# xrange和range类似，但是是返回一个迭代器，性能更优。但是在python3中，xrange和range是合并了的
# xrange(start, stop[, step])
# 搅乱训练集，让其排序顺序发生变化
temp = list(zip(x_train, y_train))
random.shuffle(temp)
x_train[:], y_train[:] = zip(*temp)
# 按照小样本数量划分训练集
mini_batches_x = [
x_train[k:k+mini_batch_size]
for k in range(0, n_train, mini_batch_size)
]
mini_batches_y = [
y_train[k:k+mini_batch_size]
for k in range(0, n_train, mini_batch_size)
]
# print("mini_batches_x=", np.array(mini_batches_x).shape)
# print("mini_batches_y= ", np.array(mini_batches_y).shape)
for (mini_x,mini_y) in zip(mini_batches_x,mini_batches_y):
# 根据每个小样本来更新w和b
self.updata_mini_batch(mini_x,mini_y,eta)
# 输出测试每轮结束后，神经网络的准确度
if True:
print("Epoch {0}: {1} / {2}".format(j, self.evaluate(x_test,y_test), n_test))
#pdb.set_trace()
else:
# 若没有数据集的情况下，则表示该次训练结束？
print("Epoch {0} complete".format(j))
def updata_mini_batch(self, mini_x, mini_y, eta):
# 作用：根据每个小样本来更新w和b
# param mini_batch：一部分的样本
# parma eta：学习率
# 根据b和w的大小分别创建全为0的矩阵
nabla_b = [np.zeros(b.shape) for b in self.biases]
nabla_w = [np.zeros(w.shape) for w in self.weights]
for (x, y) in zip(mini_x, mini_y):
# 根据样本中的每一个输入x的其输出y，计算w和b的偏导数
delta_nabla_b, delta_nabla_w = self.backprop(x,y)
# 累加存储偏导值
nabla_b = [nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
nabla_w = [nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
# 根据累加的偏导值更新w和b，这里因为使用了小样本，所以eta要除以小样本的长度
self.weights = [w-(eta/mini_x.shape[0])*nw
for w, nw in zip(self.weights, nabla_w)]
self.biases = [b-(eta/mini_x.shape[0])*nb
for b, nb in zip(self.biases, nabla_b)]
def backprop(self, x, y):
# 作用：计算偏导数
# param x：输入一个样本
# param y：样本对应的值
# return：b，w
##print("x= ", x.shape)
##print("y= ",y.shape)
nabla_b = [np.zeros(b.shape) for b in self.biases]
nabla_w = [np.zeros(w.shape) for w in self.weights]
# 前向传输
activation = x.reshape(x.shape[0], 1) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
# 存储每层的神经元的值的矩阵
activations = [x]
# 存储每个未经过sigmoid计算的神经元的值
zs = []
for b,w in zip(self.biases, self.weights):
z = np.dot(w, activation) + b
zs.append(z)
activation = sigmoid(z)
activations.append(activation)
# 求 δ 的值
delta = self.cost_derivative(activations[-1], y) * sigmoid_prime(zs[-1])
nabla_b[-1] = delta
# 乘于前一层的输出值
nabla_w[-1] = np.dot(delta, activations[-2].transpose())
for l in range(2, self.num_layers):
# 从倒数第l层开始更新
# 下面这里利用 l+1 层的 δ 值来计算 l 的 δ 值
z = zs[-l]
sp = sigmoid_prime(z)
delta = np.dot(self.weights[-l+1].transpose(), delta) * sp
nabla_b[-l] = delta
nabla_w[-l] = np.dot(delta, activations[-l-1].reshape(1, activations[-l-1].shape[0])) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
return (nabla_b, nabla_w)
def evaluate(self, x_test, y_test):
# 获得预测结果
test_results = [(np.argmax(self.feedforward(x)), y)
for (x, y) in zip(x_test, y_test)]
# 返回正确识别的个数
# print(test_results)
res = sum(int(x == y) for (x, y) in test_results)
print("This is the result: ", res)
return res
def cost_derivative(self, output_activations, y):
# 作用：二次损失函数
# :param output_activations:
# :param y: 一个数字，如6
# :return:
e = np.zeros((10,1)) # 矩阵间的减法，需要对齐。在相应数字的位置置1，其余是0
e[y] = 1.0
return (output_activations-e)
其中：
def __init__(self, sizes):
# 作用：初始化w和b
# param sizes: list 类型，存储每层神经网络的神经元数目
# 例如sizes=[2,3,2]表示输出层有两个神经元，隐藏层有3个神经元，输出层有2个神经元
# num_layers代表神经网络数目
self.num_layers = len(sizes)
self.sizes = sizes
# 初始输入层，随机产生每层中的y个神经元的biase值（0,1）
self.biases = [np.random.randn(y, 1) for y in sizes[1:]]
self.weights = [np.random.randn(y,x) for x, y in zip(sizes[:-1],sizes[1:])]
# biases和weights是共用的参数。只有一套
# biases和weights其实是一个不规则的二维数组。如果输入是3层，那么数组有2行
# 第一行30个维度，第二行10个维度
# 这个weights的初始方法算是学到了。假设sizes是[a1,a2,a3,a4]
# 那么输出是（a1,a2）(a2,a3) (a3,a4)
这是初始化函数，对我们输入的神经网络格式，我们初始化相应的数据结构。主要是对权重w和偏置b进行初始化。
def feedforward(self, x):
# 作用：向前传输计算每个神经元的值
# param a：输入值
# return：计算后每个神经元的值
a = x.reshape(x.shape[0],1) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
for b, w in zip(self.biases, self.weights):
# 加权求和再加上biase
a = sigmoid(np.dot(w,a)+b)
# dot是矩阵乘法的意思
return a
这部分是向前输入的步骤，对应以下的公式：
def SGD(self, x_train, y_train, x_test, y_test, epochs, mini_batch_size, eta):
# 作用：随机梯度下降
# param training_data：输入的训练集
# param epochs：迭代次数
# param mini_batch_size：小样本数量
# param eta：学习率
# param test_data：测试数据集
n_test = x_test.shape[0]
n_train = x_train.shape[0] # n_train是总的训练样本数目
for j in range(epochs):
# xrange和range类似，但是是返回一个迭代器，性能更优。但是在python3中，xrange和range是合并了的
# xrange(start, stop[, step])
# 搅乱训练集，让其排序顺序发生变化
temp = list(zip(x_train, y_train))
random.shuffle(temp)
x_train[:], y_train[:] = zip(*temp)
# 按照小样本数量划分训练集
mini_batches_x = [
x_train[k:k+mini_batch_size]
for k in range(0, n_train, mini_batch_size)
]
mini_batches_y = [
y_train[k:k+mini_batch_size]
for k in range(0, n_train, mini_batch_size)
]
# print("mini_batches_x=", np.array(mini_batches_x).shape)
# print("mini_batches_y= ", np.array(mini_batches_y).shape)
for (mini_x,mini_y) in zip(mini_batches_x,mini_batches_y):
# 根据每个小样本来更新w和b
self.updata_mini_batch(mini_x,mini_y,eta)
# 输出测试每轮结束后，神经网络的准确度
if True:
print("Epoch {0}: {1} / {2}".format(j, self.evaluate(x_test,y_test), n_test))
#pdb.set_trace()
else:
# 若没有数据集的情况下，则表示该次训练结束？
print("Epoch {0} complete".format(j))
这部分其实是对数据集进行“拆分”，拆分成一小份一小份，我们每次根据这一份份数据来对权重w和偏置b进行更新。
def updata_mini_batch(self, mini_x, mini_y, eta):
# 作用：根据每个小样本来更新w和b
# param mini_batch：一部分的样本
# parma eta：学习率
# 根据b和w的大小分别创建全为0的矩阵
nabla_b = [np.zeros(b.shape) for b in self.biases]
nabla_w = [np.zeros(w.shape) for w in self.weights]
for (x, y) in zip(mini_x, mini_y):
# 根据样本中的每一个输入x的其输出y，计算w和b的偏导数
delta_nabla_b, delta_nabla_w = self.backprop(x,y)
# 累加存储偏导值
nabla_b = [nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
nabla_w = [nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
# 根据累加的偏导值更新w和b，这里因为使用了小样本，所以eta要除以小样本的长度
self.weights = [w-(eta/mini_x.shape[0])*nw
for w, nw in zip(self.weights, nabla_w)]
self.biases = [b-(eta/mini_x.shape[0])*nb
for b, nb in zip(self.biases, nabla_b)]
这部分的内容就是用小样本来更新参数。对应下面的式子：
def backprop(self, x, y):
# 作用：计算偏导数
# param x：输入一个样本
# param y：样本对应的值
# return：b，w
##print("x= ", x.shape)
##print("y= ",y.shape)
nabla_b = [np.zeros(b.shape) for b in self.biases]
nabla_w = [np.zeros(w.shape) for w in self.weights]
# 前向传输
activation = x.reshape(x.shape[0], 1) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
# 存储每层的神经元的值的矩阵
activations = [x]
# 存储每个未经过sigmoid计算的神经元的值
zs = []
for b,w in zip(self.biases, self.weights):
z = np.dot(w, activation) + b
zs.append(z)
activation = sigmoid(z)
activations.append(activation)
# 求 δ 的值
delta = self.cost_derivative(activations[-1], y) * sigmoid_prime(zs[-1])
nabla_b[-1] = delta
# 乘于前一层的输出值
nabla_w[-1] = np.dot(delta, activations[-2].transpose())
for l in range(2, self.num_layers):
# 从倒数第l层开始更新
# 下面这里利用 l+1 层的 δ 值来计算 l 的 δ 值
z = zs[-l]
sp = sigmoid_prime(z)
delta = np.dot(self.weights[-l+1].transpose(), delta) * sp
nabla_b[-l] = delta
nabla_w[-l] = np.dot(delta, activations[-l-1].reshape(1, activations[-l-1].shape[0])) # 将数据由(x,)转换成(x,1)。否则矩阵计算会出错
return (nabla_b, nabla_w)
这部分是计算偏导，除了上面两个式子，还有下面的式子：
def evaluate(self, x_test, y_test):
# 获得预测结果
test_results = [(np.argmax(self.feedforward(x)), y)
for (x, y) in zip(x_test, y_test)]
# 返回正确识别的个数
# print(test_results)
res = sum(int(x == y) for (x, y) in test_results)
print("This is the result: ", res)
return res
def cost_derivative(self, output_activations, y):
# 作用：二次损失函数
# :param output_activations:
# :param y: 一个数字，如6
# :return:
e = np.zeros((10,1)) # 矩阵间的减法，需要对齐。在相应数字的位置置1，其余是0
e[y] = 1.0
return (output_activations-e)
这两个函数就比较简单了，分别是计算预测结果和二次损失函数。
到这里类的方法就结束了，下面是通用的方法。
def sigmoid(z):
#作用：求 sigmoid 函数的值
# param z:
# return:
return 1.0/(1.0+np.exp(-z))
def sigmoid_prime(z):
# 作用：求 sigmoid 函数的导数
# :param z:
# :return:
return sigmoid(z)*(1-sigmoid(z))
这是关于sigmoid函数的实现和其导数。
# %pdb on
if __name__ == "__main__":
(x_train, y_train), (x_test, y_test) = mnist.load_data() # 加载数据集，原本的数据格式是60000*28*28。三维的形式
x_train = x_train.reshape(x_train.shape[0], x_train.shape[1]*x_train.shape[2])
x_test = x_test.reshape(x_test.shape[0], x_test.shape[1]*x_test.shape[2])
y_train = y_train.reshape(y_train.shape[0],1)
y_test = y_test.reshape(y_test.shape[0],1)
x_train = x_train.astype('float32') # 将数据类型转为float32，因为很多时候我们用numpy从文本文件读取数据作为numpy的数组，默认的dtype是float64
x_test = x_test.astype('float32')
x_train /= 255
x_test /= 255
print(x_train.shape) # (60000, 784)
print(y_train.shape) # (60000, 1)
print(x_test.shape) # (10000, 784)
print(y_test.shape) # (10000, 1)
print("load data finish")
pdb.set_trace() # 设置中断，输入c继续运行
net = Network([784, 30, 10])
print("net.w1= ", np.array(net.biases[0]).shape)
print("net.w2= ", np.array(net.biases[1]).shape)
print("net.w1= ", np.array(net.weights[0]).shape)
print("net.w2= ", np.array(net.weights[1]).shape)
net.SGD(x_train, y_train, x_test, y_test, 30, 10, 3.0)
这部分是主函数。主要是加载数据集，并进行一定的转换。然后声明网络，并进行训练和分类。结果大致如下：
(60000, 784)
(60000, 1)
(10000, 784)
(10000, 1)
load data finish
> <ipython-input-10-84adde9c4133>(23)<module>()
-> net = Network([784, 30, 10])
(Pdb) c
net.w1= (30, 1)
net.w2= (10, 1)
net.w1= (30, 784)
net.w2= (10, 30)
This is the result: 7983
Epoch 0: 7983 / 10000
This is the result: 9144
Epoch 1: 9144 / 10000
This is the result: 9276
Epoch 2: 9276 / 10000
This is the result: 9296
Epoch 3: 9296 / 10000
This is the result: 9285
Epoch 4: 9285 / 10000
This is the result: 9300
Epoch 5: 9300 / 10000
This is the result: 9296
Epoch 6: 9296 / 10000
This is the result: 9303
Epoch 7: 9303 / 10000
This is the result: 9289
Epoch 8: 9289 / 10000
This is the result: 9291
Epoch 9: 9291 / 10000
This is the result: 9278
Epoch 10: 9278 / 10000
This is the result: 9289
Epoch 11: 9289 / 10000
This is the result: 9288
Epoch 12: 9288 / 10000
This is the result: 9278
Epoch 13: 9278 / 10000
This is the result: 9273
Epoch 14: 9273 / 10000
This is the result: 9275
Epoch 15: 9275 / 10000
This is the result: 9277
Epoch 16: 9277 / 10000
This is the result: 9264
Epoch 17: 9264 / 10000
This is the result: 9267
Epoch 18: 9267 / 10000
This is the result: 9264
Epoch 19: 9264 / 10000
This is the result: 9267
Epoch 20: 9267 / 10000
This is the result: 9260
Epoch 21: 9260 / 10000
This is the result: 9254
Epoch 22: 9254 / 10000
This is the result: 9253
Epoch 23: 9253 / 10000
This is the result: 9258
Epoch 24: 9258 / 10000
This is the result: 9262
Epoch 25: 9262 / 10000
This is the result: 9260
Epoch 26: 9260 / 10000
This is the result: 9260
Epoch 27: 9260 / 10000
This is the result: 9263
Epoch 28: 9263 / 10000
This is the result: 9257
Epoch 29: 9257 / 10000
小结：以上是使用BP算法实现MNIST数据集的代码啦，具体的可以参考我的github。其中需要比较注意的矩阵的相加和相乘。区别(x,1) 和 (x,)之间的区别，否则会导致矩阵运算错误。并且BP算法的思想比较简单，建议大家弄懂推导的过程。
如有不足之处，还望大家多多指教。
扫一扫
6368
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
