深度学习 梯度下降 反向传播 数字识别

深度学习实战-从原理到实践

Posted by TheMatrix on 2020-01-28

现当今机器学习/深度学习技术在某些具体垂直领域已被大量广泛应用到现实世界中，已经不再像前几年那么“火热”，与之对应的各类深度学习框架也是“百花齐放,百家争鸣”，框架终究只是个工具，不过简化了从“零”开始复杂繁琐的工作，让很多普通人都可以快速入门。本博客不单纯完成一个任务，也不涉及过多理论推导，而是真正体会到算法工作一步步原理，逐步实现，岂不乐乎？

本篇博客以经典的MNIST手写数字识别为例，逐步一步步实现通用的深度学习网络模型架构，不调用任何第三方库和框架，使用matlab进行快速搭建、训练和测试。程序中所涉及的理论知识及使用的变量名严格按照DNN神经网络的反向更新（BP）、卷积神经网络(CNN)反向传播算法 这两篇博客的符号和公式进行。MNIST手写数字包含60000张训练图片，10000张测试图片，图片大小为28×28，灰度图像，官网给出的是四个二进制存储的文件,分别为训练和测试的数据集和标签文件。假设读者已经明白所给链接博客的理论知识（不清楚可以参考更多文后的文献和程序代码中给的链接），我们接下来进行下面的具体实现。

网络架构设计

考虑到网络简单和易用性，根据MNIST数据集特点，设计了四层网络层，分别为conv+relu+meanPool、conv+relu+dropout、conv+relu+dropout、conv+softmax四个连续模块层。每个模块层进行串连连接，不涉及跳层结构，其中第一个模块层的conv是CNN卷积操作，卷积核大小为9×9×1×20,四维矩阵（h×w×c_in×c_out），即20个卷积核；第二个及以后模块的conv其实是全连接层，因为前面的层每个神经元都会与下一层所有神经元进行全连接，神经元个数分别人工设定为95、45、10。这三层是用BP反向传播原理进行参数更新，第一层是按照CNN的反向传播算法进行更新。网络示意图如下,示意图中只画了CNN层和一个最后模块层：
img

网络训练工作流程

1、数字图像训练集预处理：具体是60000张图像和标签通过二进制方式逐步读入，然后归一化为28×28×1×numsImg大小，[0,1]范围，float类型数据；标签为numsImg×1 大小,0~9数字，float类型数据；
2、超参数设定：具体是网络权重学习速率alpha=0.01，动量因子beta = 0.01, dropout层丢弃因子ratio = 0.01，总循环代数epoch=2,批处理大小batchSize = 10;权重和梯度初始化，注意大小和范围；这些超参数主要根据经验确定。
3、反向传播算法+梯度下降算法，迭代更新参数寻优。
训练主文件为“train.m”,
训练预览部分数字图像情况：
img1
训练迭代过程：
img2

网络测试工作流程

待网络训练完毕，测试流程与训练流程大致相似，也需要进行正向传播，但也有些区别，测试过程不涉及到网络参数更新，训练的一些因素也需要冻结屏蔽。
1、数字图像测试集预处理：具体是10000张图像和标签通过二进制方式逐步读入，然后归一化为28×28×1×numsImg大小，[0,1]范围，float类型数据；标签为numsImg×1 大小,0~9数字，float类型数据；
2、超参数学习速率，动量因子，dropout等需要屏蔽，屏蔽方式见"Predict.m"内容。
3、计算统计测试整体准确率。
测试主文件为“test.m”,
img3
测试集总体准确率：
img4

code

** talk is cheap,show me the code!**
项目代码地址
train.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
%% 1、数据MNISTData预处理
imageFileNameTrain = 'G:\MNIST\train-images.idx3-ubyte';
labelFileNameTtrain = 'G:\MNIST\train-labels.idx1-ubyte';
nClasses = 10;
order = 0:9;

[X_Train,Label_Train] = processMNISTdata(imageFileNameTrain,labelFileNameTtrain);
Label_Train_h = onehot(Label_Train,nClasses,order);

montage(X_Train(:,:,:,1:9))
title(['Ground Truth:',num2str(Label_Train(1:9)')]);

%% 2、超参数设定
alpha=0.01;
beta =0.01;
ratio = 0.01;
epoch = 2;
batchSize = 10;

W1=1e-2*randn(9,9,1,20);% 第一层为CNN卷积，9*9大小，输入通道为1，输出通道设定为20
W2=(2*rand(95,2000)-1)/20;% 第二层为BP卷积层，设定神经元个数为95，2000为第一层特征融合计算得出的
W3=(2*rand(45,95)-1)/10; % 第三层为BP卷积层，设定神经元个数为45
W4=(2*rand(10,45)-1)/5; % 第四层为BP卷积层，设定神经元个数为10，因为直接与10个数字对应
mmt1 = zeros(size(W1));% 带动量的W1梯度，稳定训练
mmt2 = zeros(size(W2));
mmt3 = zeros(size(W3));
mmt4 = zeros(size(W4));

%% 3、反向传播算法+梯度下降算法，迭代更新参数寻优
numCorrect = 0;% 累计预测正确的样本个数
numAll = 0;% 累计输入网络的样本个数
accuracy = 0;
for i = 1:epoch
    [~, ~, ~,numsImgs] = size(X_Train);
    for    idx= 1:batchSize:numsImgs
        %% batch data
        batchInds = idx:min(idx+batchSize-1,numsImgs);
        batchX = X_Train(:,:,:,batchInds);% 28*28*batchSize
        batchY = Label_Train_h(:,batchInds);% 10*batchSize
        [~,~,~,bs] = size(batchX);
        % show images and true labels
        % montage(batchX);
        % title(['Ground Truth:',num2str(Label_Train(batchInds)')]);
        
        %% 第一层：CNN卷积+relu+池化层
        Z1= ConvLayer(batchX,W1);
        Y1=ReLULayer(Z1);
        A1 =PoolLayer(Y1);
        
        %% 第二层：BP全连接层 ＋relu+ DropoutLayer
        y1=reshape(A1,[],batchSize);
        Z2=W2*y1;
        Y2=ReLULayer(Z2);
        [A2,mask2] = DropoutLayer(Y2, ratio);
        
        %% 第三层：BP全连接层 ＋ relu+ DropoutLayer
        Z3=W3*A2;
        Y3=ReLULayer(Z3);
        [A3,mask3] =  DropoutLayer(Y3, ratio);
        
        %% 第四层：BP全连接层+softmax
        Z4 = W4*A3;
        P = SoftmaxLayer(Z4); % 10*batchSize 大小
        % 计算训练集当前准确率
        [~,ind] = max(P);
        predict_L = onehot(ind-1,nClasses,order);
        for idx_img = 1:bs
            isEqual = predict_L(:,idx_img)==batchY(:,idx_img);
            numCorrect = numCorrect+all(isEqual);
        end
        numAll = numAll+bs;
        accuracy = numCorrect/numAll;
        fprintf('第%d epoch，第%d/%d代总体训练集准确率为：%.2f\n',i,floor(idx/batchSize),floor(numsImgs/batchSize),accuracy);
        
        %% 递推求误差
        % 求第四层误差，即最后一层误差
        e4 = batchY-P;
        
        % 求第三层误差
        delta4 = e4;
        e3=W4'*delta4;
        delta3=mask3.*e3;% DropoutLayer求导 https://blog.csdn.net/oBrightLamp/article/details/84105097
        delta3 = (Z3>0).*delta3;% Relu求导
        % 求第二层误差
        e2=W3'*delta3;
        delta2=mask2.*e2;
        delta2=(Z2>0).*delta2;
        % 求第一层误差
        e1=W2'*delta2;
        e1 = reshape(e1,size(A1));
        avg_e = e1/4; % avg pool层求导参考 https://blog.csdn.net/qq_21190081/article/details/72871704
        E1 = repelem(avg_e,2,2);
        delta1=(Z1>0).*E1;
        
        % 求取梯度并更新W1,按照CNN卷积层的误差反向传播原理https://www.cnblogs.com/pinard/p/6494810.html
        [~,~,~,numFilters]=size(W1);
        for idx_f =1:numFilters
            dW1 = zeros(size(W1));% 9*9*1*20
            for idx_img = 1:bs
                dW1(:,:,:,idx_img)=alpha* convn(batchX(:,:,:,idx_img),rot90(delta1(:,:,idx_f,idx_img),2),'valid');
            end
            dW1 = mean(dW1,4); % batch 梯度方向均值
            S = size(W1);S(end) = 1;
            dW1= reshape(dW1,S);% 保持维度
            mmt1(:,:,:,idx_f)= dW1 + beta*mmt1(:,:,:,idx_f);
            W1(:,:,:,idx_f)=W1(:,:,:,idx_f) + mmt1(:,:,:,idx_f); % 带动量形式的mini-batch梯度下降算法
        end
        % 求取梯度并更新W2,按照BP反向传播原理
        dW2=alpha*delta2*y1';% 因为第一层是CNN，卷积后的特征是全连接，故y1是A1的变换形式
        mmt2 = dW2 + beta*mmt2;
        W2   = W2 + mmt2;
        % 求取梯度并更新W3,按照BP反向传播原理
        dW3=alpha*delta3*A2';
        mmt3 = dW3 + beta*mmt3;
        W3   = W3 + mmt3;
        % 求取梯度并更新W4,按照BP反向传播原理
        dW4=alpha*delta4*A3';
        mmt4 = dW4 + beta*mmt4;
        W4   = W4 + mmt4;
        
    end
    % 每个epoch完后进行保存
    save(['model_epoch',num2str(i),'.mat'], 'W1','W2','W3','W4');  
end

test.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
%% 测试集测试
load model_epoch2.mat
imageFileNameTest = 'G:\MNIST\t10k-images.idx3-ubyte';
labelFileNameTest = 'G:\MNIST\t10k-labels.idx1-ubyte';
nClasses = 10;
order = 0:9;
% 预处理
[X_Test,Label_Test] = processMNISTdata(imageFileNameTest,labelFileNameTest);
Label_true = onehot(Label_Test,nClasses,order);% 10*numImgs
% 随意查看前16个数字情况
montage(X_Test(:,:,:,1:16))
title(['Ground Truth:',num2str(Label_Test(1:16)')]);

%% 所有样本整体准确度
[~,~,~,numImgs] = size(X_Test);
predict_L=Predict(W1,W2,W3,W4,X_Test);
numCorrect = 0;
for idx_img = 1:numImgs
    isEqual = predict_L(:,idx_img)==Label_true(:,idx_img);
    numCorrect = numCorrect+all(isEqual);
end
numAll = numImgs;
acc = numCorrect/numAll;

fprintf('测试集整体准确率为：%.5f\n',acc);

ConvLayer.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function Z = ConvLayer(imgs, W)
% 功能：对图像imgs进行卷积操作,W为卷积核
% 输入：imgs 待卷积的图像或特征，存储顺序为 img_height*img_width*img_channel*batchSize,
%            四维数组，float,[0,1]范围;
%      W 为卷积核,大小为 kernel_height*hernel_width*kernel_channel*out_filters,
%         四维数组，float类型,kernel_channel必须与img_channel相等才可以卷积!
% 输出：
%      Z为卷积后的特征图，维度大小为(img_height-kernel_height+1)*(img_width-hernl_width+1)*out_filters*batchSize
%        四维数组，float类型;
% 注意：此为图像卷积操作，边界无填充，支持多维卷积操作
%
%  author:cuixingxing 2020.1.25
% email:cuixingxing150@gmail.com
%

[img_height, img_width, img_channel,batchSize] = size(imgs);
[kernel_height, hernl_width, kernel_channel, out_filters] = size(W);
assert(kernel_channel==img_channel);

Z = zeros(img_height-kernel_height+1,img_width-hernl_width+1,out_filters,batchSize);
for i = 1:batchSize
    for j = 1:out_filters
        W = rot90(W,2);% 用convn函数做卷积需要把卷积核旋转180度,只转第一维，第二维度
        Z(:,:,j,i) = convn(imgs(:,:,:,i),W(:,:,:,j),'valid');
    end
end
end

DropoutLayer.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function [out_features,mask] = DropoutLayer(in_features, ratio)
% 功能：以ratio的比例丢弃神经元
% 输入：in_features,输入特征
%       ratio 丢弃比例，[0,1]之间
% 输出：out_features，输出特征,大小与in_features一致
%       mask，随机掩码，大小与in_features一致
% 参考：https://zhuanlan.zhihu.com/p/38200980
%      https://blog.csdn.net/oBrightLamp/article/details/84105097
%
% author:cuixingxing 2020.1.27
% email:cuixingxing150@gmail.com
%

all_nums = numel(in_features);
drop_nums = round(all_nums*ratio);
idxs = randperm(all_nums,drop_nums);
mask = ones(size(in_features));
mask(idxs) = 0;

out_features = in_features.*mask./(1-ratio);
end

onehot.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function onehotMatrix = onehot(labels,nClasses,order)
% 产生独热编码矩阵，根据order顺序返回onehot矩阵，其中每列为一个one hot标签
% 输入：labels 包含所有类别标签的类别，1*nums或者nums*1单个类别或者字符向量
%      nClasses 包含类别总数
%      order 长度为numLabels的标签顺序，向量，labels中元素来自于order
% 输出：onehotMatrix
%       独热编码矩阵，nClasses*nums大小，每列中有且仅有一个1，其余为0
% example:
%  labels = [0,3,2,8];
%  nClasses = 10;
%  order = [0,1,2,3,4,5,6,7,8,9];
%  onehotMatrix = onehot(labels,nClasses,order);
%  onehotMatrix = [1  0  0  0
%                  0  0  0  0
%                  0  0  1  0
%                  0  1  0  0
%                  0  0  0  0
%                  0  0  0  0
%                  0  0  0  0
%                  0  0  0  0
%                  0  0  0  1
%                  0  0  0  0]
%
% author:cuixingxing 2020.1.27
% cuixingxing150@gmail.com
%

assert(nClasses==length(order));
assert(all(ismember(labels,order)));
labels = categorical(labels);
order = categorical(order);

E = eye(nClasses);
nums = numel(labels);
indexs = [];
for i = 1:nums
    ind = find(labels(i)==order);
    indexs = [indexs;ind];
end
onehotMatrix = E(:,indexs);

PoolLayer.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function out = PoolLayer(in)
% 功能：2*2 平均池化层,stride = 2,pad = 0
% 输入：in 特征层 in_height*in_width*in_channel*batchSize,
%            四维数组，float;
% 输出：out 特征层 out_height/2*out_width/2*out_channel*batchSize,
%            四维数组，float;
%
% 注意：out_channel == in_channel
%
% author:cuixingxing 2020.1.26
% email:cuixingxing150@gmail.com
%
[in_height,in_width,in_channel,batchSize] = size(in);
temp = zeros(in_height-2+1,in_width-2+1,in_channel,batchSize);
W = ones(2,2,in_channel)/4;
for i = 1:batchSize
    for j = 1:in_channel
       temp(:,:,j,i) = convn(in(:,:,j,i),W(:,:,j),'valid');
    end
end
out = temp(1:2:end,1:2:end,:,:);
end

Predict.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function predict_L = Predict(W1,W2,W3,W4,X_Test)
% 功能：计算测试集准确率
% 输入：W1,W2,W3,W4为网络权重矩阵
%      X_Test 测试数据，存储顺序为 img_height*img_width*img_channel*numImgs,
%            四维数组，float类型,[0,1]范围;
% 输出：predict_L 预测标签，onehot标签，每列为一个标签
%
%  author:cuixingxing 2020.1.27
% email:cuixingxing150@gmail.com
%

%% 预处理
X = X_Test; % 28*28*1*numImgs， float类型，[0,1]范围
[~,~,~,numImgs] = size(X);

%% forward
%% 第一层：CNN卷积+ReLULayer+池化层
Z1= ConvLayer(X,W1);
Y1=ReLULayer(Z1);
A1 =PoolLayer(Y1);

%% 第二层：全连接层 ＋ReLULayer+ dropout
y1=reshape(A1,[],numImgs);
Z2=W2*y1;
Y2=ReLULayer(Z2);
% [A2,~] = DropoutLayer(Y2, 0.01);

%% 第三层：全连接层 ＋ ReLULayer+ dropout
Z3=W3*Y2;
Y3=ReLULayer(Z3);
% [A3,~] =  DropoutLayer(Y3, 0.01);

%% 第四层：全连接层+softmax
Z4 = W4*Y3;
% 计算训练集准确率
P = SoftmaxLayer(Z4); % 10*batchSize 大小
[~,ind] = max(P);
nClasses = 10;
order = 0:9;
predict_L = onehot(ind-1,nClasses,order);

end

preProcess.m

1
2
3
4
5
6
7
8
9
function out = preProcess(in)
% 功能：预处理输入图像
% 输入：in 为H*W*numImgs 大小，uint8类型图像
% 输出：out 为H*W*C*numImgs大小，float类型[0,1]范围图像,mnist数字图像，其中C=1
%
[H,W,numImgs] = size(in);
out = im2single(in);
out = reshape(out,[H,W,1,numImgs]);
end

processMNISTdata.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function [Imgs,Labels] = processMNISTdata(imageFileName,labelFileName)
% PROCESSMNISTDATA功能：解析官方提供的二进制数字识别MNIST数据集
%  输入：imageFileName 解压后的图像文件名
%       labelFileName 解压后的标签文件名
% 输出：Imgs 为numRows*numCols*1*numImages 大小的float类型图像数据，[0,1]范围
%       Labels为numImages*1大小的double标签数据
% 参考：https://ww2.mathworks.cn/help/stats/visualize-high-dimensional-data-using-t-sne.html
%
%  author:cuixingxing 2020.1.25
% email:cuixingxing150@gmail.com
%
[fileID,errmsg] = fopen(imageFileName,'r','b');
if fileID < 0
    error(errmsg);
end
%%
% First read the magic number. This number is 2051 for image data, and
% 2049 for label data
magicNum = fread(fileID,1,'int32',0,'b');
assert(magicNum == 2051);

%%
% Then read the number of images, number of rows, and number of columns
numImages = fread(fileID,1,'int32',0,'b');
numRows = fread(fileID,1,'int32',0,'b');
numCols = fread(fileID,1,'int32',0,'b');

%%
% Read the image data
Imgs = fread(fileID,inf,'unsigned char=>uint8');
%%
% Reshape the data to array X
Imgs = reshape(Imgs,numCols,numRows,numImages);
Imgs = permute(Imgs,[2 1 3]);
Imgs = preProcess(Imgs); % 28*28*1*numImgs， float类型，[0,1]范围

%%
% Close the file
fclose(fileID);
%%
% Similarly, read the label data.
[fileID,errmsg] = fopen(labelFileName,'r','b');
if fileID < 0
    error(errmsg);
end
magicNum = fread(fileID,1,'int32',0,'b');
assert(magicNum == 2049);
numItems = fread(fileID,1,'int32',0,'b');

Labels = fread(fileID,inf,'unsigned char=>double');
Labels = reshape(Labels,numItems,1);% numItem*1 
fclose(fileID);

ReLULayer.m

1
2
3
function y = ReLULayer(x)
  y = max(0, x);
end

SoftmaxLayer.m

1
2
3
4
5
6
7
8
function y = SoftmaxLayer(x)
%% 功能：softmax转概率
% 输入：x 10*numsSamples大小矩阵
% 输出：y 10*numsSamples概率矩阵，每列为一个样本概率
%
  ex = exp(x);
  y  = ex ./ sum(ex,1);
end

Reference

神经网络BP反向传播算法原理和详细推导流程
梯度下降(Gradient Descent)简析及matlab实现
深度学习笔记（3）——CNN中一些特殊环节的反向传播

FEATURED TAGS
深度学习 梯度下降 反向传播 数字识别
FRIENDS