引言
这篇文章的内容是导师给的小作业，其实是在很久之前就已经完成了的，现在有些忘记生疏了，故写一篇文章总结一下。
这篇文章中，将具体阐述如何通过计算机用简单的神经网络算法来实现手写数字的识别，并简单阐述背后的基本原理。文章的后面可以看到， 经过训练后的神经网络对手写数字的识别可以达到94%~95的准确率。
参考资料
《机器学习》 周志华著
Neural Networks and Deep Learning
关于这个题目其实已经有更详细的文章来讲述，如果你愿意建议花时间看一下neuralnetworksanddeeplearning.com/chap1.html。我主要是以自己的理解来写下这篇文章。文章的代码部分是用MATLAB来实现的。文章结构是由“论述+代码”组成的，每一个小标题的文字论述后面都会有具体的代码实现。
1.我们的目标
首先需要明确一下我们的目标：识别手写数字。这里我们采用的是来自MINST的手写数据集，由60000个数字图片以及对应的人工识别数字组成。我们拿一些出来看一下：
我们很容易知道这些数字图片的规格为28x28像素。
如何下载并在MATLAB中导入这些数据呢？
这里为了方便，我直接给出下载地址：识别标签 图片数据
每一张图片都是灰度图片，即每一个像素点都是由0-255来描述它的黑色浓度的（姑且这么叫吧）。我们将他们下载好了之后需要将导入MATLAB,并做一些简单的处理。包括将28x28的矩阵转换为784x1的矩阵，归一化处理将每一个像素的值转换到0~1之间（很容易想到除以255就可以了），将数据文件中的一些描述性文字去除。
这里就不再赘述具体实现，直接调用下面的代码就可以将数据导入。
这是导入标签数据集的function代码：
function labels = loadMNISTLabels(filename)
%loadMNISTLabels returns a [number of MNIST images]x1 matrix containing
%the labels for the MNIST images
fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);
magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2049, ['Bad magic number in ', filename, '']);
numLabels = fread(fp, 1, 'int32', 0, 'ieee-be');
labels = fread(fp, inf, 'unsigned char');
assert(size(labels,1) == numLabels, 'Mismatch in label count');
fclose(fp);
end
导入图片数据集的function代码：
function images = loadMNISTImages(filename)
%loadMNISTImages returns a 28x28x[number of MNIST images] matrix containing
%the raw MNIST images
fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);
magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2051, ['Bad magic number in ', filename, '']);
numImages = fread(fp, 1, 'int32', 0, 'ieee-be');
numRows = fread(fp, 1, 'int32', 0, 'ieee-be');
numCols = fread(fp, 1, 'int32', 0, 'ieee-be');
images = fread(fp, inf, 'unsigned char');
images = reshape(images, numCols, numRows, numImages);
images = permute(images,[2 1 3]);
fclose(fp);
% Reshape to #pixels x #examples
images = reshape(images, size(images, 1) * size(images, 2), size(images, 3));
% Convert to double and rescale to [0,1]
images = double(images) / 255;
end
2.如何进行识别
我们已经知道对于，识别这样的数字对人脑来说十分简单。但怎么让计算机也可以对这样的手写数字进行识别呢？
正如标题所说这里我们用到的是神经网络的算法，具体而言是一个单隐层前馈网络。
接下来我们具体地来论述这个所谓的“单隐层前馈网络”是什么意思，它代表了什么。以及对这个网络结构为什么能给出一张手写数字的识别结果建立一些直觉上的认识。
在生物神经网络中，每个神经元与其他的神经元相连接，当一个神经元的电位大于一定的阈值时，它就会被激活从而“兴奋”，并且向其他神经元发送化学物质，引起其他神经元电位状态的改变。
现在我们先抛开上面那张乱七八糟的图片不管，先来看一个简单的东西，它的名字叫做“感知器”（Perceptron），是一种人造神经元。
一个感知器可以有多个输入，一个输出。上图是一个拥有3个输入一个输出的感知器。它的工作原理是这样的：先把我们的每一个输入乘以一个权重，将他们累加起来，当累加的值超过了阈值（threshold）时，感知器的输出为1，否则为0。用数学语言描述如下：
如何直观地来感受这个感知器呢？举个例子吧，假设我们需要用一个感知器来判断一部手机是不是你想要的？我们为这个感知器设计3个输入参数：X1是屏幕大小，X2是屏占比，X3是手机重量，输出1时为你想要的，0时则不是。
我们可以通过调整权重和阈值来进行决策，也就是判断一部手机是否是你想要的。比如你比较在乎的是屏幕和重量，而对屏占比不怎么看重。那么这就意味着屏幕和和重量的权重会比较大，而屏占比的权重比较小。也就是说，当我们给定了这一个感知器的所有权重和阈值时，拿一部手机到你面前，你是否想要这一部手机的决策可以由这个感知器给出。
现在比如我们想要用单个感知器识别一张图片是否为9。
通过观察上面这张图片,我们人工地给这样一个感知器大概的描述：
我们给图片中的黑色部分的位置一个较高的权值，说明这些位置的黑色浓度对一张图片是不是9的决策影响很大；给白色部分的权值赋一个很小的值，说明当这些位置的黑色浓度对一张图片是不是9的决策几乎没有影响。这就是感知器做出决策的过程，而我们上面这个“人工观察”的过程可以由计算机学习多张图片来完成，具体过程不再赘述。
需要注意：感知器只有输出层神经元进行激活函数处理，即只拥有一层功能神经元，其学习能力非常有限。显然单个感知器的是无法完全地解决我们的目标：“识别手写数字”的。这里我们采用一种称为“多层前馈神经网络”的结构，如下所示：
每层神经元与下一层神经元全互连，神经元之间不存在同层连接，也不存在跨层连接，每一层的输出作为下一层的输入，这样的神经网络结构称之为“多层前馈神经网络”。第一层为输入层，最后一层为输出层，中间的层称之为隐层。显然，对比我们之前的单感知器结构，只是多了中间的隐层而已。而我们的功能神经元便是由隐层和输出层的神经元来组成的。
需要注意的是：输入层是不对数据进行处理的，它只是把所有的数据原封不动地引入而已。
现在回到我们的手写识别问题，我们的每一张图片是由784个像素组成的，所以我们的一张图片应该对应了784输入层神经元。而我们的识别范围是0-9，因此需要完全表征识别的数字至少需要10个输出层神经元。现在还剩下隐层的神经元没有确定。实际上隐层的层数和神经元的数目是不固定的，且隐层层数越多，神经元数目越大识别效果越好，但也就意味着学习时间会加长。
这里我们采用的是单隐层，15个隐层神经元的网络结构。如开头所述的，这是一个单隐层前馈网络。
而实际上，除了网络结构更复杂外，有一些细节也必须强调：
1.为了更加清晰和直观地来描述每一个神经元的数学模型，我们用bias = - threshold来替代我们上面描述的阈值。
2.我们将神经元的输出定义为总输出值（乘以了权重后的输入值累加）加上bias后，再经过一个激活函数后所输出的值。
当f为阶跃函数即sgn()函数时，这个神经元就是我们上面所说的感知器。而一般我们采用的激活函数f为sigmoid()函数，相较于前者它具有光滑连续的特点。由sigmoid函数定义的神经元也称为S神经元。
两个函数的定义：
现在我们新的网络结构已经定义完成，我们试着来直观地理解一下当我们的网络训练完成后它是怎么给出识别结果的。通过观察网络结构发现其实我们的10个输出神经元的输出决策过程是独立的，我们输出的这样一个10x1的向量其实是对输入的判断是否是“0”“1”...“9”的十种识别结果（或者说可能性）合在一起罢了。以“5”对应的输出神经元为例，当我们的784个像素具有这样的一种特征的时候（比如输入神经元的第5，8，17....687的权重比较大）它是“5”的可能性就很大，对应了我们的输出层的第6个神经元的输出很大（接近1）。这里有一个问题是我们并无法直观上理解中间的隐层到底做了什么，所以我们上面所谓的直观理解其实是不准确的，它不像没有隐层的感知器那样很容易被直观理解。这里我们暂时认为中间的隐层找到了“5”对应的像素特征从而导出了正确的结果。
至此，我们已经用文字完整地论述了神经网络的结构。接下来便是代码的实现。
神经网络的代码实现：
观察上面那张神经网络的结构你可能会想到用类的方式来实现，即定义一个神经元类，然后产生一个神经元矩阵。在matlab中我们更一般的做法是用矩阵的方式，这和上面的类的方式不太一样，但也是很好理解的。
我们的神经网络结构由输入层，隐层，输出层组成。输入层只是把我们的数据原封不动的搬入而已，不需要任何参数。隐层需要30个bias，以及30X784个连接输入层和隐层的weights；同理，输出层有10个bias，10X30个连接输出层和隐层的weights。当我们给出上面所述的所有参数后，其实我们的神经网络就已经完全确定了！我们用的是randn函数进行初始化。
初始化网络中的各种参数：
%使用randn初始化biases和weights
biases_h=randn(30,1);%隐层biases
biases_o=randn(10,1);%输出层biases
weights_ih=randn(30,784);%连接输入层和隐层的权值
weights_ho=randn(10,30);%连接隐层和输出层的权值
然后我们定义我们的output函数，它可以导出一层神经元在接受到上一层神经元的输入后的输出结果，它的定义是分两步进行的：
1.定义sigmoid函数：
function y = sigmoid( a )
%SIGMOID sigmoid函数定义
y=1./(1+exp(-a));
end
2.定义output函数:
function y = output( i,w,b )
%OUTPUT 返回某一层神经元的输出
% 参数顺序为i,w,b;i代表input，w代表weights，b代表biases
temp=w*i+b;
y=sigmoid(temp);
end
在完成上面几步之后，我们定义神经网络的工作就全部完成了（只用了10行代码）。比如我们现在想要得到一张784X1的图像输入经过神经网络后的输出结果，只需要调用两次output函数就可以了。
举例：如何得到一张784X1图像的输出：
Input=zeros(784,1);%输入数据，假设是元素全为0的矩阵
Output_hidden=output(Input,weights_ih,biases_h);%得到隐层的输出
Output=output(Output_hidden,weights_ho,biases_o);%得到输出层的输出，也就是神经网络的输出
3.如何训练我们的神经网络
神经网络的训练过程其实就是根据训练数据来调整神经元之间的连接权(weights)，以及每一个功能神经元的阈值(biases)；我们学习到的内容蕴涵在连接权和阈值中。
我们这里采用的是一种称为“误差逆传播”的算法（简称BP算法），由于前馈神经网络多采用BP算法来训练，故前馈神经网络也称为BP网络。我们接下来看一下它是怎么推导出来的以及具体如何运用它来训练我们的网络的。（后面是一些数学推导，如果你很讨厌的话可以直接跳过）
这里由于简书上对公式的支持不太友好，我们接下来的论述直接用图片+文字给出：（这里如果看不太清楚可以查看原图）
咳咳，看图片真是累人。我也是，搞了半天才把这几张图片给完成。别看有好几张图，其实涉及的数学很简单，无非就是求偏导数，链式法则。你只要会求偏导数就可以独立推导出BP算法的表达式出来。好了，如果你看懂了上面所论述的，那现在训练网络的方法你就已经完全了解了。
在展示代码之前有两点点必须强调：
1.我们这里采用的是随机批梯度下降。也就是每一次先随机打乱我们的训练数据集，然后每次取10个出来计算更新量，接着更新所有的参数。10个训练案例的累计更新量等于单独每个训练案例的更新量累加再除以10（这里关于累计更新量的定义是有依据的）。
2.我们得到的数据集是60000张图片和对应的标签。这里我们取50000个作为训练集，剩下10000个作为测试集。
BP算法的代码实现
这里代码是分几步来进行的：
1.首先将MINST数据集分成训练集和测试集两个部分。其中训练集和测试集都是由图片矩阵和标签所组成的。
2.随机打乱我们的训练集，按顺序每次取10个训练案例。
3.由这10个训练案例计算出各参数的更新量，然后更新所有的参数。
4.循环2-3直至遍历完训练集的每一个训练案例。
5.循环2-4直至准确率足够高。
这里为了便于观察训练过程中我们识别准确率的变化，我们在步骤4之后做一些小小的改变。
5.遍历完一次训练集后，利用测试集来验证我们网络的准确率，输出“正确识别个数/10000”。
6.循环2-5步骤30次。
代码里还有一些细节没有提到，如果你足够仔细的话，配合注释理解起来应该没有问题。只要找准矩阵位置和所对应的参数就好了。
将MINST数据集分成训练集和测试集
%获取数据并将数据拆分成训练数据和测试数据
%前50000个为训练数据，后10000个为测试数据
IMGS=loadMNISTImages('train-images.idx3-ubyte');
TRAIN_IMGS=IMGS(:,1:50000);
TEST_IMGS=IMGS(:,50001:60000);
LABS=loadMNISTLabels('train-labels.idx1-ubyte');
TRAIN_LABS=LABS(1:50000,:);
TEST_LABS=LABS(50001:60000,:);
%%将识别标签转换为向量形式
%比如7则转换为[0,0,0,0,0,0,0,1,0,0]
%TRAIN_LABS
TEMP=zeros(5003000,10);
for i=1:50000
TEMP(i,TRAIN_LABS(i,1)+1)=1;
end
TRAIN_LABS=TEMP.';
%TEST_LABS
TEMP=zeros(10000,10);
for i=1:10000
TEMP(i,TEST_LABS(i,1)+1)=1;
end
TEST_LABS=TEMP.';
利用BP算法训练神经网络并输出准确率
%这里定义了out_o是给后面测试集的验证使用的，先分配了空间
out_o=zeros(10,10000);
%%现在我们的网络已经初始化完成，下面采用mini批BP算法来学习
%在这里我们设置每一批的数目为10
step=3;%设置学习率
for j=1:30 %j代表训练次数，暂设置最大训练数为100
%%下面的FOR循环为对数据集的一次mini批BP算法应用
%首先随机打乱我们的数据集
r=randperm(50000);
TRAIN_IMGS=TRAIN_IMGS(:,r);
TRAIN_LABS=TRAIN_LABS(:,r);
for i=0:4999
%计算一个mini-batch的输出层神经元的梯度项G，即把该mini-batch的每一项求得的G相加
G=zeros(10,1);%分配空间
E=zeros(30,1);%分配空间
%给各个参数的delta值分配空间，delta即变化量
weights_ho_delta=zeros(10,30);
biases_o_delta=zeros(10,1);
weights_ih_delta=zeros(30,784);
biases_h_delta=zeros(30,1);
for l=(i*10+1):((i+1)*10)
%首先由一个样本获得输出
tem_i=TRAIN_IMGS(:,l);
out_h= output(tem_i,weights_ih,biases_h);
tem_o=output(out_h,weights_ho,biases_o);
tem_labs=TRAIN_LABS(:,l);
%计算输出层神经元的梯度项G
G=tem_o.*(1-tem_o).*(tem_labs-tem_o);
%计算隐层神经元的梯度项E
E=out_h.*(1-out_h).*((weights_ho.')*G);
%获得各个参数的delta值，原理是把mini-batch里的每一样本求delta值然后累加，再/10
weights_ho_delta=weights_ho_delta+(step/10)*(G*(out_h.'));
biases_o_delta=biases_o_delta-(step/10)*G;
weights_ih_delta=weights_ih_delta+(step/10)*(E*(tem_i.'));
biases_h_delta=biases_h_delta-(step/10)*E;
end
%%接下来由上面求出的E和G更新参数
%更新隐层到输出层的权值
weights_ho=weights_ho+weights_ho_delta;
%更新输出层的biases
biases_o=biases_o+biases_o_delta;
%更新输入层和隐层的权值
weights_ih=weights_ih+weights_ih_delta;
%更新隐层的biases
biases_h=biases_h+biases_h_delta;
end
%%每一次迭代运用测试数据集来验证准确率，输出测试集的正确输出数目
%对于神经网络的输出做如下处理：
%由测试集的输入通过神经网络获得输出（10*10000的矩阵）,取每一列的最大值为1，其他为0。
%比如某一列为[0.1,0.3,0.2,0.2,0.5,0.7,0.8,0.1,0,0].'将会被转成[0,0,0,0,0,0,1,0,0,0].'
for i=1:10000
%获得一个样本输入的神经网络输出
tem_i=TEST_IMGS(:,i);
out_h= output(tem_i,weights_ih,biases_h);
tem_o=output(out_h,weights_ho,biases_o);
out_o(:,i)=tem_o;
%转换（即将最大值所在的位置取1，其他取0）
tem_row=1;
tem_max=out_o(tem_row,i);
for k=1:10
if out_o(k,i)>tem_max
tem_row=k;
tem_max=out_o(tem_row,i);
end
end
for k=1:10
if k==tem_row
out_o(k,i)=1;
else
out_o(k,i)=0;
end
end
end
%与标准输出对比，获得准确率
sum=0;
for i=1:10000
if out_o(:,i)==TEST_LABS(:,i)
sum=sum+1;
end
end
fprintf('迭代%d: %d/10000\n',j,sum);
end
4.完整的代码以及一些程序运行结果
loadMNISTLabels.m
function labels = loadMNISTLabels(filename)
%loadMNISTLabels returns a [number of MNIST images]x1 matrix containing
%the labels for the MNIST images
fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);
magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2049, ['Bad magic number in ', filename, '']);
numLabels = fread(fp, 1, 'int32', 0, 'ieee-be');
labels = fread(fp, inf, 'unsigned char');
assert(size(labels,1) == numLabels, 'Mismatch in label count');
fclose(fp);
end
loadMNISTImages.m
function images = loadMNISTImages(filename)
%loadMNISTImages returns a 28x28x[number of MNIST images] matrix containing
%the raw MNIST images
fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);
magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2051, ['Bad magic number in ', filename, '']);
numImages = fread(fp, 1, 'int32', 0, 'ieee-be');
numRows = fread(fp, 1, 'int32', 0, 'ieee-be');
numCols = fread(fp, 1, 'int32', 0, 'ieee-be');
images = fread(fp, inf, 'unsigned char');
images = reshape(images, numCols, numRows, numImages);
images = permute(images,[2 1 3]);
fclose(fp);
% Reshape to #pixels x #examples
images = reshape(images, size(images, 1) * size(images, 2), size(images, 3));
% Convert to double and rescale to [0,1]
images = double(images) / 255;
end
sigmoid.m
function y = sigmoid( a )
%SIGMOID sigmoid函数定义
y=1./(1+exp(-a));
end
output.m
function y = output( i,w,b )
%OUTPUT 返回某一层神经元的输出
% 参数顺序为i,w,b;i代表input，w代表weights，b代表biases
temp=w*i+b;
y=sigmoid(temp);
end
main.m
%获取数据并将数据拆分成训练数据和测试数据
%前50000个为训练数据，后10000个为测试数据
IMGS=loadMNISTImages('train-images.idx3-ubyte');
TRAIN_IMGS=IMGS(:,1:50000);
TEST_IMGS=IMGS(:,50001:60000);
LABS=loadMNISTLabels('train-labels.idx1-ubyte');
TRAIN_LABS=LABS(1:50000,:);
TEST_LABS=LABS(50001:60000,:);
%%将识别标签转换为向量形式
%比如7则转换为[0,0,0,0,0,0,0,1,0,0]
%TRAIN_LABS
TEMP=zeros(5003000,10);
for i=1:50000
TEMP(i,TRAIN_LABS(i,1)+1)=1;
end
TRAIN_LABS=TEMP.';
%TEST_LABS
TEMP=zeros(10000,10);
for i=1:10000
TEMP(i,TEST_LABS(i,1)+1)=1;
end
TEST_LABS=TEMP.';
%定义网络大小w2
NetSize=[784,30,10];
%使用randn初始化biases和weights
biases_h=randn(30,1);%隐层biases
biases_o=randn(10,1);%输出层biases
weights_ih=randn(30,784);%连接输入层和隐层的权值
weights_ho=randn(10,30);%连接隐层和输出层的权值
%这里定义了out_o是给后面测试集的验证使用的，先分配了空间
out_o=zeros(10,10000);
%%现在我们的网络已经初始化完成，下面采用mini批BP算法来学习
%在这里我们设置每一批的数目为10
step=3;%设置学习率
for j=1:30 %j代表训练次数，暂设置最大训练数为100
%%下面的FOR循环为对数据集的一次mini批BP算法应用
%首先随机打乱我们的数据集
r=randperm(50000);
TRAIN_IMGS=TRAIN_IMGS(:,r);
TRAIN_LABS=TRAIN_LABS(:,r);
for i=0:4999
%计算一个mini-batch的输出层神经元的梯度项G，即把该mini-batch的每一项求得的G相加
G=zeros(10,1);%分配空间
E=zeros(30,1);%分配空间
%给各个参数的delta值分配空间，delta即变化量
weights_ho_delta=zeros(10,30);
biases_o_delta=zeros(10,1);
weights_ih_delta=zeros(30,784);
biases_h_delta=zeros(30,1);
for l=(i*10+1):((i+1)*10)
%首先由一个样本获得输出
tem_i=TRAIN_IMGS(:,l);
out_h= output(tem_i,weights_ih,biases_h);
tem_o=output(out_h,weights_ho,biases_o);
tem_labs=TRAIN_LABS(:,l);
%计算输出层神经元的梯度项G
G=tem_o.*(1-tem_o).*(tem_labs-tem_o);
%计算隐层神经元的梯度项E
E=out_h.*(1-out_h).*((weights_ho.')*G);
%获得各个参数的delta值，原理是把mini-batch里的每一样本求delta值然后累加，再/10
weights_ho_delta=weights_ho_delta+(step/10)*(G*(out_h.'));
biases_o_delta=biases_o_delta-(step/10)*G;
weights_ih_delta=weights_ih_delta+(step/10)*(E*(tem_i.'));
biases_h_delta=biases_h_delta-(step/10)*E;
end
%%接下来由上面求出的E和G更新参数
%更新隐层到输出层的权值
weights_ho=weights_ho+weights_ho_delta;
%更新输出层的biases
biases_o=biases_o+biases_o_delta;
%更新输入层和隐层的权值
weights_ih=weights_ih+weights_ih_delta;
%更新隐层的biases
biases_h=biases_h+biases_h_delta;
end
%%每一次迭代运用测试数据集来验证准确率，输出测试集的正确输出数目
%对于神经网络的输出做如下处理：
%由测试集的输入通过神经网络获得输出（10*10000的矩阵）,取每一列的最大值为1，其他为0。
%比如某一列为[0.1,0.3,0.2,0.2,0.5,0.7,0.8,0.1,0,0].'将会被转成[0,0,0,0,0,0,1,0,0,0].'
for i=1:10000
%获得一个样本输入的神经网络输出
tem_i=TEST_IMGS(:,i);
out_h= output(tem_i,weights_ih,biases_h);
tem_o=output(out_h,weights_ho,biases_o);
out_o(:,i)=tem_o;
%转换（即将最大值所在的位置取1，其他取0）
tem_row=1;
tem_max=out_o(tem_row,i);
for k=1:10
if out_o(k,i)>tem_max
tem_row=k;
tem_max=out_o(tem_row,i);
end
end
for k=1:10
if k==tem_row
out_o(k,i)=1;
else
out_o(k,i)=0;
end
end
end
%与标准输出对比，获得准确率
sum=0;
for i=1:10000
if out_o(:,i)==TEST_LABS(:,i)
sum=sum+1;
end
end
fprintf('迭代%d: %d/10000\n',j,sum);
end
输出结果图
可以看到，我们的网络在训练之后可以达到接近95%的准确率，这个结果已经相当令人满意了对吧？毕竟我们只用了不到100行代码（注释不算）就完成了我们的所有工作。
如果你自己手写两个数字进去识别一般情况下也是可以得出正确结果的，这里代码也很简单，就不再赘述了。
好了，那这篇文章到这里就结束了，如果看完你能对神经网络的结构有一些更深入的了解那我会很高兴的。之后有时间会更新一下一些更好理解的可视化动态图，敬请期待吧。
在这里可以直接下载我的程序以及数据集http://pan.baidu.com/s/1jI9WyJs。
