Deep Learning模型之：CNN的反向求导及练习

[1]Deep learning简介

[2]Deep Learning训练过程

[3]Deep Learning模型之：CNN卷积神经网络推导和实现

[4]Deep Learning模型之：CNN的反向求导及练习

[5]Deep Learning模型之：CNN卷积神经网络（一）深度解析CNN

[6]Deep Learning模型之：CNN卷积神经网络（二）文字识别系统LeNet-5

[7]Deep Learning模型之：CNN卷积神经网络（三）CNN常见问题总结

前言：

CNN作为DL中最成功的模型之一，有必要对其更进一步研究它。虽然在前面的博文Stacked CNN简单介绍中有大概介绍过CNN的使用，不过那是有个前提的：CNN中的参数必须已提前学习好。而本文的主要目的是介绍CNN参数在使用bp算法时该怎么训练，毕竟CNN中有卷积层和下采样层，虽然和MLP的bp算法本质上相同，但形式上还是有些区别的，很显然在完成CNN反向传播前了解bp算法是必须的。本文的实验部分是参考斯坦福UFLDL新教程UFLDL：Exercise: Convolutional Neural Network里面的内容。完成的是对MNIST

 

　　实验基础：

　　CNN反向传播求导时的具体过程可以参考论文Notes on Convolutional Neural Networks, Jake Bouvrie，该论文讲得很全面，比如它考虑了pooling层也加入了权值、偏置值及非线性激发(因为这2种值也需要learn)，对该论文的解读可参考zouxy09的博文CNN卷积神经网络推导和实现。除了bp算法外，本人认为理解了下面4个子问题，基本上就可以弄懂CNN的求导了（bp算法这里就不多做介绍，网上资料实在是太多了），另外为了讲清楚一些细节过程，本文中举的例子都是简化了一些条件的，且linux下文本和公式编辑太难弄了，文中难免有些地方会出错，大家原谅下。

　　首先我们来看看CNN系统的目标函数，设有样本(xi, yi)共m个，CNN网络共有L层，中间包含若干个卷积层和pooling层，最后一层的输出为f(xi)，则系统的loss表达式为(对权值进行了惩罚，一般分类都采用交叉熵形式)：

 　　

 

　　问题一：求输出层的误差敏感项。

　　现在只考虑个一个输入样本(x, y)的情形，loss函数和上面的公式类似是用交叉熵来表示的，暂时不考虑权值规则项，样本标签采用one-hot编码，CNN网络的最后一层采用softmax全连接(多分类时输出层一般用softmax)，样本(x,y)经过CNN网络后的最终的输出用f(x)表示，则对应该样本的loss值为:

 　　

　　其中f(x)的下标c的含义见公式：

 　　

　　因为x通过CNN后得到的输出f(x)是一个向量，该向量的元素值都是概率值，分别代表着x被分到各个类中的概率，而f(x)中下标c的意思就是输出向量中取对应c那个类的概率值。

　　采用上面的符号，可以求得此时loss值对输出层的误差敏感性表达式为：

 　　

　　其中e(y)表示的是样本x标签值的one-hot表示，其中只有一个元素为1,其它都为0.

　　其推导过程如下（先求出对输出层某个节点c的误差敏感性，参考Larochelle关于DL的课件:Output layer gradient）,求出输出层中c节点的误差敏感值：

 　　

　　由上面公式可知，如果输出层采用sfotmax，且loss用交叉熵形式，则最后一层的误差敏感值就等于CNN网络输出值f(x)减样本标签值e(y),即f(x)-e(y),其形式非常简单，这个公式是不是很眼熟？很多情况下如果model采用MSE的loss，即loss=1/2*(e(y)-f(x))^2，那么loss对最终的输出f(x)求导时其结果就是f(x)-e(y),虽然和上面的结果一样，但是大家不要搞混淆了，这2个含义是不同的，一个是对输出层节点输入值的导数(softmax激发函数)，一个是对输出层节点输出值的导数(任意激发函数）。而在使用MSE的loss表达式时，输出层的误差敏感项为(f(x)-e(y)).*f(x)’，两者只相差一个因子。

　　这样就可以求出第L层的权值W的偏导数：

 　　

　　输出层偏置的偏导数：

 　　

　　上面2个公式的e(y)和f(x)是一个矩阵，已经把所有m个训练样本考虑进去了,每一列代表一个样本。

 

　　问题二：当接在卷积层的下一层为pooling层时，求卷积层的误差敏感项。

　　假设第l(小写的l，不要看成数字’1’了)层为卷积层，第l+1层为pooling层，且pooling层的误差敏感项为：   ,卷积层的误差敏感项为：  , 则两者的关系表达式为：

 　　

　　这里符号●表示的是矩阵的点积操作，即对应元素的乘积。卷积层和unsample()后的pooling层节点是一一对应的，所以下标都是用j表示。后面的符号表示的是第l层第j个节点处激发函数的导数(对节点输入的导数)。

　　其中的函数unsample()为上采样过程，其具体的操作得看是采用的什么pooling方法了。但unsample的大概思想为：pooling层的每个节点是由卷积层中多个节点(一般为一个矩形区域)共同计算得到，所以pooling层每个节点的误差敏感值也是由卷积层中多个节点的误差敏感值共同产生的，只需满足两层见各自的误差敏感值相等，下面以mean-pooling和max-pooling为例来说明。

　　假设卷积层的矩形大小为4×4, pooling区域大小为2×2, 很容易知道pooling后得到的矩形大小也为2*2（本文默认pooling过程是没有重叠的，卷积过程是每次移动一个像素，即是有重叠的，后续不再声明）,如果此时pooling后的矩形误差敏感值如下：

 　　

　　则按照mean-pooling，首先得到的卷积层应该是4×4大小，其值分布为(等值复制)：

 　　

　　因为得满足反向传播时各层见误差敏感总和不变，所以卷积层对应每个值需要平摊（除以pooling区域大小即可，这里pooling层大小为2×2=4)），最后的卷积层值

分布为：

 　　

　　mean-pooling时的unsample操作可以使用matalb中的函数kron()来实现，因为是采用的矩阵Kronecker乘积。C=kron(A, B)表示的是矩阵B分别与矩阵A中每个元素相乘，然后将相乘的结果放在C中对应的位置。比如：

　　　　

　　如果是max-pooling，则需要记录前向传播过程中pooling区域中最大值的位置，这里假设pooling层值1,3,2,4对应的pooling区域位置分别为右下、右上、左上、左下。则此时对应卷积层误差敏感值分布为：

　　 

　　当然了，上面2种结果还需要点乘卷积层激发函数对应位置的导数值了，这里省略掉。

 

　　问题三：当接在pooling层的下一层为卷积层时，求该pooling层的误差敏感项。

　　假设第l层(pooling层)有N个通道，即有N张特征图，第l+1层(卷积层)有M个特征，l层中每个通道图都对应有自己的误差敏感值，其计算依据为第l+1层所有特征核的贡献之和。下面是第l+1层中第j个核对第l层第i个通道的误差敏感值计算方法：

　　

　　符号★表示的是矩阵的卷积操作，这是真正意义上的离散卷积，不同于卷积层前向传播时的相关操作，因为严格意义上来讲，卷积神经网络中的卷积操作本质是一个相关操作，并不是卷积操作，只不过它可以用卷积的方法去实现才这样叫。而求第i个通道的误差敏感项时需要将l+1层的所有核都计算一遍，然后求和。另外因为这里默认pooling层是线性激发函数，所以后面没有乘相应节点的导数。

　　举个简单的例子，假设拿出第l层某个通道图，大小为3×3，第l+1层有2个特征核，核大小为2×2，则在前向传播卷积时第l+1层会有2个大小为2×2的卷积图。如果2个特征核分别为：

 　　    

　　反向传播求误差敏感项时，假设已经知道第l+1层2个卷积图的误差敏感值：

 　　 

　　离散卷积函数conv2()的实现相关子操作时需先将核旋转180度(即左右翻转后上下翻转)，但这里实现的是严格意义上的卷积，所以在用conv2()时，对应的参数核不需要翻转（在有些toolbox里面，求这个问题时用了旋转，那是因为它们已经把所有的卷积核都旋转过，这样在前向传播时的相关操作就不用旋转了。并不矛盾）。且这时候该函数需要采用’full’模式，所以最终得到的矩阵大小为3×3,（其中3=2+2-1）,刚好符第l层通道图的大小。采用’full’模式需先将第l+1层2个卷积图扩充，周围填0,padding后如下：

 　　       

　　扩充后的矩阵和对应的核进行卷积的结果如下情况：

 　　

　　 

　　可以通过手动去验证上面的结果，因为是离散卷积操作，而离散卷积等价于将核旋转后再进行相关操作。而第l层那个通道的误差敏感项为上面2者的和，呼应问题三，最终答案为：

 　　

　　那么这样问题3这样解的依据是什么呢？其实很简单，本质上还是bp算法，即第l层的误差敏感值等于第l+1层的误差敏感值乘以两者之间的权值，只不过这里由于是用了卷积，且是有重叠的，l层中某个点会对l+1层中的多个点有影响。比如说最终的结果矩阵中最中间那个0.3是怎么来的呢？在用2×2的核对3×3的输入矩阵进行卷积时，一共进行了4次移动，而3×3矩阵最中间那个值在4次移动中均对输出结果有影响，且4次的影响分别在右下角、左下角、右上角、左上角。所以它的值为2×0.2+1×0.1+1×0.1-1×0.3=0.3, 建议大家用笔去算一下，那样就可以明白为什么这里可以采用带’full’类型的conv2()实现。

 

　　问题四：求与卷积层相连那层的权值、偏置值导数。

　　前面3个问题分别求得了输出层的误差敏感值、从pooling层推断出卷积层的误差敏感值、从卷积层推断出pooling层的误差敏感值。下面需要利用这些误差敏感值模型中参数的导数。这里没有考虑pooling层的非线性激发，因此pooling层前面是没有权值的，也就没有所谓的权值的导数了。现在将主要精力放在卷积层前面权值的求导上(也就是问题四)。

　　假设现在需要求第l层的第i个通道，与第l+1层的第j个通道之间的权值和偏置的导数，则计算公式如下：

 　　

　　其中符号⊙表示矩阵的相关操作，可以采用conv2()函数实现。在使用该函数时，需将第l+1层第j个误差敏感值翻转。

　　例如，如果第l层某个通道矩阵i大小为4×4,如下：

 　　

　　第l+1层第j个特征的误差敏感值矩阵大小为3×3,如下：

 　　

　　很明显，这时候的特征Kij导数的大小为2×2的，且其结果为：

 　　

　　而此时偏置值bj的导数为1.2 ，将j区域的误差敏感值相加即可(0.8+0.1-0.6+0.3+0.5+0.7-0.4-0.2=1.2)，因为b对j中的每个节点都有贡献，按照多项式的求导规则(和的导数等于导数的和)很容易得到。

　　为什么采用矩阵的相关操作就可以实现这个功能呢？由bp算法可知，l层i和l+1层j之间的权值等于l+1层j处误差敏感值乘以l层i处的输入，而j中某个节点因为是由i+1中一个区域与权值卷积后所得，所以j处该节点的误差敏感值对权值中所有元素都有贡献，由此可见，将j中每个元素对权值的贡献(尺寸和核大小相同)相加，就得到了权值的偏导数了(这个例子的结果是由9个2×2大小的矩阵之和)，同样，如果大家动笔去推算一下，就会明白这时候为什么可以用带’valid’的conv2()完成此功能。

 

　　实验总结：

1. 卷积层过后，可以先跟pooling层，再通过非线性传播函数。也可以是先通过非线性传播函数再经过pooling层。
2. CNN的结构本身就是一种规则项。
3. 实际上每个权值的学习率不同时优化会更好。
4. 发现Ng以前的ufldl中教程里面softmax并没有包含偏置值参数，至少他给的start code里面没有包含，严格来说是错误的。
5. 当输入样本为多个时，bp算法中的误差敏感性也是一个矩阵。每一个样本都对应有自己每层的误差敏感性。
6. 血的教训啊，以后循环变量名不能与终止名太相似，一不小心引用下标是就弄错，比如for filterNum = 1:numFilters 时一不小心把下标用numFilters去代替了，花了大半天去查这个debug.

　　7.  matlab中conv2()函数在卷积过程中默认是每次移动一个像素，即重叠度最高的卷积。

 

　　实验结果：

　　按照网页教程UFLDL：Exercise: Convolutional Neural Network和UFLDL：Optimization: Stochastic Gradient Descent，对MNIST数据库进行识别，完成练习中YOU CODE HERE部分后，该CNN网络的识别率为：

　　95.76%

　　只采用了一个卷积层+一个pooling层+一个softmax层。卷积层的特征个数为20,卷积核大小为9×9, pooling区域大小为2×2.

　　在进行实验前，需下载好本实验的标准代码：https://github.com/amaas/stanford_dl_ex。

　　然后在common文件夹放入下载好的MNIST数据库，见http://yann.lecun.com/exdb/mnist/.注意MNIST文件名需要和代码中的保持一致。

 

　　实验代码：

cnnTrain.m: 

[plain] view plain copy

1. </pre><pre code_snippet_id="456994" snippet_file_name="blog_20140823_2_375615" name="code" class="plain" style="font-size: 13px;">%% Convolution Neural Network Exercise  
2.   
3. %  Instructions  
4. %  ------------  
5. %   
6. %  This file contains code that helps you get started in building a single.  
7. %  layer convolutional nerual network. In this exercise, you will only  
8. %  need to modify cnnCost.m and cnnminFuncSGD.m. You will not need to   
9. %  modify this file.  
10.   
11. %%======================================================================  
12. %% STEP 0: Initialize Parameters and Load Data  
13. %  Here we initialize some parameters used for the exercise.  
14.   
15. % Configuration  
16. imageDim = 28;  
17. numClasses = 10;  % Number of classes (MNIST images fall into 10 classes)  
18. filterDim = 9;    % Filter size for conv layer,9*9  
19. numFilters = 20;   % Number of filters for conv layer  
20. poolDim = 2;      % Pooling dimension, (should divide imageDim-filterDim+1)  
21.   
22. % Load MNIST Train  
23. addpath ./common/;   
24. images = loadMNISTImages('./common/train-images-idx3-ubyte');  
25. images = reshape(images,imageDim,imageDim,[]);  
26. labels = loadMNISTLabels('./common/train-labels-idx1-ubyte');  
27. labels(labels==0) = 10; % Remap 0 to 10  
28.   
29. % Initialize Parameters,theta=(2165,1),2165=9*9*20+20+100*20*10+10  
30. theta = cnnInitParams(imageDim,filterDim,numFilters,poolDim,numClasses);  
31.   
32. %%======================================================================  
33. %% STEP 1: Implement convNet Objective  
34. %  Implement the function cnnCost.m.  
35.   
36. %%======================================================================  
37. %% STEP 2: Gradient Check  
38. %  Use the file computeNumericalGradient.m to check the gradient  
39. %  calculation for your cnnCost.m function.  You may need to add the  
40. %  appropriate path or copy the file to this directory.  
41.   
42. DEBUG=false;  % set this to true to check gradient  
43. %DEBUG = true;  
44. if DEBUG  
45.     % To speed up gradient checking, we will use a reduced network and  
46.     % a debugging data set  
47.     db_numFilters = 2;  
48.     db_filterDim = 9;  
49.     db_poolDim = 5;  
50.     db_images = images(:,:,1:10);  
51.     db_labels = labels(1:10);  
52.     db_theta = cnnInitParams(imageDim,db_filterDim,db_numFilters,...  
53.                 db_poolDim,numClasses);  
54.       
55.     [cost grad] = cnnCost(db_theta,db_images,db_labels,numClasses,...  
56.                                 db_filterDim,db_numFilters,db_poolDim);  
57.       
58.   
59.     % Check gradients  
60.     numGrad = computeNumericalGradient( @(x) cnnCost(x,db_images,...  
61.                                 db_labels,numClasses,db_filterDim,...  
62.                                 db_numFilters,db_poolDim), db_theta);  
63.    
64.     % Use this to visually compare the gradients side by side  
65.     disp([numGrad grad]);  
66.       
67.     diff = norm(numGrad-grad)/norm(numGrad+grad);  
68.     % Should be small. In our implementation, these values are usually   
69.     % less than 1e-9.  
70.     disp(diff);   
71.    
72.     assert(diff < 1e-9,...  
73.         'Difference too large. Check your gradient computation again');  
74.       
75. end;  
76.   
77. %%======================================================================  
78. %% STEP 3: Learn Parameters  
79. %  Implement minFuncSGD.m, then train the model.  
80.   
81. % 因为是采用的mini-batch梯度下降法，所以总共对样本的循环次数次数比标准梯度下降法要少  
82. % 很多，因为每次循环中权值已经迭代多次了  
83. options.epochs = 3;   
84. options.minibatch = 256;  
85. options.alpha = 1e-1;  
86. options.momentum = .95;  
87.   
88. opttheta = minFuncSGD(@(x,y,z) cnnCost(x,y,z,numClasses,filterDim,...  
89.                       numFilters,poolDim),theta,images,labels,options);  
90. save('theta.mat','opttheta');               
91.   
92. %%======================================================================  
93. %% STEP 4: Test  
94. %  Test the performance of the trained model using the MNIST test set. Your  
95. %  accuracy should be above 97% after 3 epochs of training  
96.   
97. testImages = loadMNISTImages('./common/t10k-images-idx3-ubyte');  
98. testImages = reshape(testImages,imageDim,imageDim,[]);  
99. testLabels = loadMNISTLabels('./common/t10k-labels-idx1-ubyte');  
100. testLabels(testLabels==0) = 10; % Remap 0 to 10  
101.   
102. [~,cost,preds]=cnnCost(opttheta,testImages,testLabels,numClasses,...  
103.                 filterDim,numFilters,poolDim,true);  
104.   
105. acc = sum(preds==testLabels)/length(preds);  
106.   
107. % Accuracy should be around 97.4% after 3 epochs  
108. fprintf('Accuracy is %f\n',acc);  

cnnConvolve.m: 

[plain] view plain copy

1. </pre><pre code_snippet_id="456994" snippet_file_name="blog_20140823_4_8183786" name="code" class="plain" style="font-size: 13px;">function convolvedFeatures = cnnConvolve(filterDim, numFilters, images, W, b)  
2. %cnnConvolve Returns the convolution of the features given by W and b with  
3. %the given images  
4. %  
5. % Parameters:  
6. %  filterDim - filter (feature) dimension  
7. %  numFilters - number of feature maps  
8. %  images - large images to convolve with, matrix in the form  
9. %           images(r, c, image number)  
10. %  W, b - W, b for features from the sparse autoencoder,传进来的w已经是矩阵的形式  
11. %         W is of shape (filterDim,filterDim,numFilters)  
12. %         b is of shape (numFilters,1)  
13. %  
14. % Returns:  
15. %  convolvedFeatures - matrix of convolved features in the form  
16. %                      convolvedFeatures(imageRow, imageCol, featureNum, imageNum)  
17.   
18. numImages = size(images, 3);  
19. imageDim = size(images, 1);  
20. convDim = imageDim - filterDim + 1;  
21.   
22. convolvedFeatures = zeros(convDim, convDim, numFilters, numImages);  
23.   
24. % Instructions:  
25. %   Convolve every filter with every image here to produce the   
26. %   (imageDim - filterDim + 1) x (imageDim - filterDim + 1) x numFeatures x numImages  
27. %   matrix convolvedFeatures, such that   
28. %   convolvedFeatures(imageRow, imageCol, featureNum, imageNum) is the  
29. %   value of the convolved featureNum feature for the imageNum image over  
30. %   the region (imageRow, imageCol) to (imageRow + filterDim - 1, imageCol + filterDim - 1)  
31. %  
32. % Expected running times:   
33. %   Convolving with 100 images should take less than 30 seconds   
34. %   Convolving with 5000 images should take around 2 minutes  
35. %   (So to save time when testing, you should convolve with less images, as  
36. %   described earlier)  
37.   
38.   
39. for imageNum = 1:numImages  
40.   for filterNum = 1:numFilters  
41.   
42.     % convolution of image with feature matrix  
43.     convolvedImage = zeros(convDim, convDim);  
44.   
45.     % Obtain the feature (filterDim x filterDim) needed during the convolution  
46.   
47.     %%% YOUR CODE HERE %%%  
48.     filter = W(:,:,filterNum);  
49.     bc = b(filterNum);  
50.       
51.     % Flip the feature matrix because of the definition of convolution, as explained later  
52.     filter = rot90(squeeze(filter),2);  
53.         
54.     % Obtain the image  
55.     im = squeeze(images(:, :, imageNum));  
56.   
57.     % Convolve "filter" with "im", adding the result to convolvedImage  
58.     % be sure to do a 'valid' convolution  
59.   
60.     %%% YOUR CODE HERE %%%  
61.     convolvedImage = conv2(im, filter, 'valid');  
62.       
63.     % Add the bias unit  
64.     % Then, apply the sigmoid function to get the hidden activation  
65.     %%% YOUR CODE HERE %%%  
66.     convolvedImage = sigmoid(convolvedImage+bc);  
67.       
68.     convolvedFeatures(:, :, filterNum, imageNum) = convolvedImage;  
69.   end  
70. end  
71.   
72. end  
73.   
74. function sigm = sigmoid(x)  
75.     sigm = 1./(1+exp(-x));  
76. end  

cnnPool.m: 

[plain] view plain copy

1. function pooledFeatures = cnnPool(poolDim, convolvedFeatures)  
2. %cnnPool Pools the given convolved features  
3. %  
4. % Parameters:  
5. %  poolDim - dimension of pooling region  
6. %  convolvedFeatures - convolved features to pool (as given by cnnConvolve)  
7. %                      convolvedFeatures(imageRow, imageCol, featureNum, imageNum)  
8. %  
9. % Returns:  
10. %  pooledFeatures - matrix of pooled features in the form  
11. %                   pooledFeatures(poolRow, poolCol, featureNum, imageNum)  
12. %       
13.   
14. numImages = size(convolvedFeatures, 4);  
15. numFilters = size(convolvedFeatures, 3);  
16. convolvedDim = size(convolvedFeatures, 1);  
17.   
18. pooledFeatures = zeros(convolvedDim / poolDim, ...  
19.         convolvedDim / poolDim, numFilters, numImages);  
20.   
21. % Instructions:  
22. %   Now pool the convolved features in regions of poolDim x poolDim,  
23. %   to obtain the   
24. %   (convolvedDim/poolDim) x (convolvedDim/poolDim) x numFeatures x numImages   
25. %   matrix pooledFeatures, such that  
26. %   pooledFeatures(poolRow, poolCol, featureNum, imageNum) is the   
27. %   value of the featureNum feature for the imageNum image pooled over the  
28. %   corresponding (poolRow, poolCol) pooling region.   
29. %     
30. %   Use mean pooling here.  
31.   
32. %%% YOUR CODE HERE %%%  
33.     % convolvedFeatures(imageRow, imageCol, featureNum, imageNum)  
34.     % pooledFeatures(poolRow, poolCol, featureNum, imageNum)  
35.     for numImage = 1:numImages  
36.         for numFeature = 1:numFilters  
37.             for poolRow = 1:convolvedDim / poolDim  
38.                 offsetRow = 1+(poolRow-1)*poolDim;  
39.                 for poolCol = 1:convolvedDim / poolDim  
40.                     offsetCol = 1+(poolCol-1)*poolDim;  
41.                     patch = convolvedFeatures(offsetRow:offsetRow+poolDim-1, ...  
42.                         offsetCol:offsetCol+poolDim-1,numFeature,numImage); %取出一个patch  
43.                     pooledFeatures(poolRow,poolCol,numFeature,numImage) = mean(patch(:));  
44.                 end  
45.             end              
46.         end  
47.     end  
48.       
49. end  

cnnCost.m:

[plain] view plain copy

1. function [cost, grad, preds] = cnnCost(theta,images,labels,numClasses,...  
2.                                 filterDim,numFilters,poolDim,pred)  
3. % Calcualte cost and gradient for a single layer convolutional  
4. % neural network followed by a softmax layer with cross entropy  
5. % objective.  
6. %                              
7. % Parameters:  
8. %  theta      -  unrolled parameter vector  
9. %  images     -  stores images in imageDim x imageDim x numImges  
10. %                array  
11. %  numClasses -  number of classes to predict  
12. %  filterDim  -  dimension of convolutional filter                              
13. %  numFilters -  number of convolutional filters  
14. %  poolDim    -  dimension of pooling area  
15. %  pred       -  boolean only forward propagate and return  
16. %                predictions  
17. %  
18. %  
19. % Returns:  
20. %  cost       -  cross entropy cost  
21. %  grad       -  gradient with respect to theta (if pred==False)  
22. %  preds      -  list of predictions for each example (if pred==True)  
23.   
24.   
25. if ~exist('pred','var')  
26.     pred = false;  
27. end;  
28.   
29. imageDim = size(images,1); % height/width of image  
30. numImages = size(images,3); % number of images  
31. lambda = 3e-3; % weight decay parameter       
32.   
33. %% Reshape parameters and setup gradient matrices  
34.   
35. % Wc is filterDim x filterDim x numFilters parameter matrix  
36. % bc is the corresponding bias  
37.   
38.   
39. % Wd is numClasses x hiddenSize parameter matrix where hiddenSize  
40. % is the number of output units from the convolutional layer  
41. % bd is corresponding bias  
42. [Wc, Wd, bc, bd] = cnnParamsToStack(theta,imageDim,filterDim,numFilters,...  
43.                         poolDim,numClasses); %the theta vector cosists wc,wd,bc,bd in order  
44.   
45. % Same sizes as Wc,Wd,bc,bd. Used to hold gradient w.r.t above params.  
46. Wc_grad = zeros(size(Wc));  
47. Wd_grad = zeros(size(Wd));  
48. bc_grad = zeros(size(bc));  
49. bd_grad = zeros(size(bd));  
50.   
51.   
52. %%======================================================================  
53. %% STEP 1a: Forward Propagation  
54. %  In this step you will forward propagate the input through the  
55. %  convolutional and subsampling (mean pooling) layers.  You will then use  
56. %  the responses from the convolution and pooling layer as the input to a  
57. %  standard softmax layer.  
58.   
59.   
60. %% Convolutional Layer  
61. %  For each image and each filter, convolve the image with the filter, add  
62. %  the bias and apply the sigmoid nonlinearity.  Then subsample the   
63. %  convolved activations with mean pooling.  Store the results of the  
64. %  convolution in activations and the results of the pooling in  
65. %  activationsPooled.  You will need to save the convolved activations for  
66. %  backpropagation.  
67. convDim = imageDim-filterDim+1; % dimension of convolved output  
68. outputDim = (convDim)/poolDim; % dimension of subsampled output  
69.   
70.   
71. % convDim x convDim x numFilters x numImages tensor for storing activations  
72. activations = zeros(convDim,convDim,numFilters,numImages);  
73.   
74.   
75. % outputDim x outputDim x numFilters x numImages tensor for storing  
76. % subsampled activations  
77. activationsPooled = zeros(outputDim,outputDim,numFilters,numImages);  
78.   
79.   
80. %%% YOUR CODE HERE %%%  
81. convolvedFeatures = cnnConvolve(filterDim, numFilters, images, Wc, bc); %前向传播,已经经过了激发函数  
82. activationsPooled = cnnPool(poolDim, convolvedFeatures);  
83.   
84. % Reshape activations into 2-d matrix, hiddenSize x numImages,  
85. % for Softmax layer  
86. activationsPooled = reshape(activationsPooled,[],numImages);  
87.   
88. %% Softmax Layer  
89. %  Forward propagate the pooled activations calculated above into a  
90. %  standard softmax layer. For your convenience we have reshaped  
91. %  activationPooled into a hiddenSize x numImages matrix.  Store the  
92. %  results in probs.  
93.   
94. % numClasses x numImages for storing probability that each image belongs to  
95. % each class.  
96. probs = zeros(numClasses,numImages);  
97.   
98. %%% YOUR CODE HERE %%%  
99. %Wd=(numClasses,hiddenSize),probs的每一列代表一个输出  
100. M = Wd*activationsPooled+repmat(bd,[1,numImages]);   
101. M = bsxfun(@minus,M,max(M,[],1));  
102. M = exp(M);  
103. probs = bsxfun(@rdivide, M, sum(M)); %why rdivide?  
104.   
105. %%======================================================================  
106. %% STEP 1b: Calculate Cost  
107. %  In this step you will use the labels given as input and the probs  
108. %  calculate above to evaluate the cross entropy objective.  Store your  
109. %  results in cost.  
110.   
111. cost = 0; % save objective into cost  
112.   
113. %%% YOUR CODE HERE %%%  
114. % cost = -1/numImages*labels(:)'*log(probs(:));  
115. % 首先需要把labels弄成one-hot编码  
116. groundTruth = full(sparse(labels, 1:numImages, 1));  
117. cost = -1./numImages*groundTruth(:)'*log(probs(:))+(lambda/2.)*(sum(Wd(:).^2)+sum(Wc(:).^2)); %加入一个惩罚项  
118. % cost = -1./numImages*groundTruth(:)'*log(probs(:));  
119.   
120.   
121. % Makes predictions given probs and returns without backproagating errors.  
122. if pred  
123.     [~,preds] = max(probs,[],1);  
124.     preds = preds';  
125.     grad = 0;  
126.     return;  
127. end;  
128.   
129. %% 将c步和d步合成在一起了  
130. %======================================================================  
131. % STEP 1c: Backpropagation  
132. %  Backpropagate errors through the softmax and convolutional/subsampling  
133. %  layers.  Store the errors for the next step to calculate the gradient.  
134. %  Backpropagating the error w.r.t the softmax layer is as usual.  To  
135. %  backpropagate through the pooling layer, you will need to upsample the  
136. %  error with respect to the pooling layer for each filter and each image.    
137. %  Use the kron function and a matrix of ones to do this upsampling   
138. %  quickly.  
139.   
140.   
141. % STEP 1d: Gradient Calculation  
142. %  After backpropagating the errors above, we can use them to calculate the  
143. %  gradient with respect to all the parameters.  The gradient w.r.t the  
144. %  softmax layer is calculated as usual.  To calculate the gradient w.r.t.  
145. %  a filter in the convolutional layer, convolve the backpropagated error  
146. %  for that filter with each image and aggregate over images.  
147.   
148.   
149. %%% YOUR CODE HERE %%%  
150. %%% YOUR CODE HERE %%%  
151. % 网络结构: images--> convolvedFeatures--> activationsPooled--> probs  
152. % Wd = (numClasses,hiddenSize)  
153. % bd = (hiddenSize,1)  
154. % Wc = (filterDim,filterDim,numFilters)  
155. % bc = (numFilters,1)  
156. % activationsPooled = zeros(outputDim,outputDim,numFilters,numImages);  
157. % convolvedFeatures = (convDim,convDim,numFilters,numImages)  
158. % images(imageDim,imageDim,numImges)  
159. delta_d = -(groundTruth-probs); % softmax layer's preactivation,每一个样本都对应有自己每层的误差敏感性。  
160. Wd_grad = (1./numImages)*delta_d*activationsPooled'+lambda*Wd;  
161. bd_grad = (1./numImages)*sum(delta_d,2); %注意这里是要求和  
162. delta_s = Wd'*delta_d; %the pooling/sample layer's preactivation  
163. delta_s = reshape(delta_s,outputDim,outputDim,numFilters,numImages);  
164.   
165.   
166. %进行unsampling操作，由于kron函数只能对二维矩阵操作，所以还得分开弄  
167. %delta_c = convolvedFeatures.*(1-convolvedFeatures).*(1./poolDim^2)*kron(delta_s, ones(poolDim));   
168. delta_c = zeros(convDim,convDim,numFilters,numImages);  
169. for i=1:numImages  
170.     for j=1:numFilters  
171.         delta_c(:,:,j,i) = (1./poolDim^2)*kron(squeeze(delta_s(:,:,j,i)), ones(poolDim));  
172.     end  
173. end  
174. delta_c = convolvedFeatures.*(1-convolvedFeatures).*delta_c;  
175.   
176.   
177. % Wc_grad = convn(images,rot90(delta_c,2,'valid'))+ lambda*Wc;  
178. for i=1:numFilters  
179.     Wc_i = zeros(filterDim,filterDim);  
180.     for j=1:numImages    
181.         Wc_i = Wc_i+conv2(squeeze(images(:,:,j)),rot90(squeeze(delta_c(:,:,i,j)),2),'valid');  
182.     end  
183.    % Wc_i = convn(images,rot180(squeeze(delta_c(:,:,i,:))),'valid');  
184.     % add penalize  
185.     Wc_grad(:,:,i) = (1./numImages)*Wc_i+lambda*Wc(:,:,i);  
186.       
187.     bc_i = delta_c(:,:,i,:);  
188.     bc_i = bc_i(:);  
189.     bc_grad(i) = sum(bc_i)/numImages;  
190. end  
191.   
192.   
193. %% Unroll gradient into grad vector for minFunc  
194. grad = [Wc_grad(:) ; Wd_grad(:) ; bc_grad(:) ; bd_grad(:)];  
195.   
196.   
197. end  
198.   
199.   
200. function X = rot180(X)  
201.     X = flipdim(flipdim(X, 1), 2);  
202. end  

minFuncSGD.m: 

[plain] view plain copy

1. function [opttheta] = minFuncSGD(funObj,theta,data,labels,...  
2.                         options)  
3. % Runs stochastic gradient descent with momentum to optimize the  
4. % parameters for the given objective.  
5. %  
6. % Parameters:  
7. %  funObj     -  function handle which accepts as input theta,  
8. %                data, labels and returns cost and gradient w.r.t  
9. %                to theta.  
10. %  theta      -  unrolled parameter vector  
11. %  data       -  stores data in m x n x numExamples tensor  
12. %  labels     -  corresponding labels in numExamples x 1 vector  
13. %  options    -  struct to store specific options for optimization  
14. %  
15. % Returns:  
16. %  opttheta   -  optimized parameter vector  
17. %  
18. % Options (* required)  
19. %  epochs*     - number of epochs through data  
20. %  alpha*      - initial learning rate  
21. %  minibatch*  - size of minibatch  
22. %  momentum    - momentum constant, defualts to 0.9  
23.   
24.   
25. %%======================================================================  
26. %% Setup  
27. assert(all(isfield(options,{'epochs','alpha','minibatch'})),...  
28.         'Some options not defined');  
29. if ~isfield(options,'momentum')  
30.     options.momentum = 0.9;  
31. end;  
32. epochs = options.epochs;  
33. alpha = options.alpha;  
34. minibatch = options.minibatch;  
35. m = length(labels); % training set size  
36. % Setup for momentum  
37. mom = 0.5;  
38. momIncrease = 20;  
39. velocity = zeros(size(theta));  
40.   
41. %%======================================================================  
42. %% SGD loop  
43. it = 0;  
44. for e = 1:epochs  
45.       
46.     % randomly permute indices of data for quick minibatch sampling  
47.     rp = randperm(m);  
48.       
49.     for s=1:minibatch:(m-minibatch+1)  
50.         it = it + 1;  
51.   
52.         % increase momentum after momIncrease iterations  
53.         if it == momIncrease  
54.             mom = options.momentum;  
55.         end;  
56.   
57.         % get next randomly selected minibatch  
58.         mb_data = data(:,:,rp(s:s+minibatch-1)); % 取出当前的mini-batch的训练样本  
59.         mb_labels = labels(rp(s:s+minibatch-1));  
60.   
61.         % evaluate the objective function on the next minibatch  
62.         [cost grad] = funObj(theta,mb_data,mb_labels);  
63.           
64.         % Instructions: Add in the weighted velocity vector to the  
65.         % gradient evaluated above scaled by the learning rate.  
66.         % Then update the current weights theta according to the  
67.         % sgd update rule  
68.           
69.         %%% YOUR CODE HERE %%%  
70.         velocity = mom*velocity+alpha*grad; % 见ufldl教程Optimization: Stochastic Gradient Descent  
71.         theta = theta-velocity;  
72.           
73.         fprintf('Epoch %d: Cost on iteration %d is %f\n',e,it,cost);  
74.     end;  
75.   
76.     % aneal learning rate by factor of two after each epoch  
77.     alpha = alpha/2.0;  
78.   
79. end;  
80.   
81. opttheta = theta;  
82.   
83. end  

computeNumericalGradient.m: 

[plain] view plain copy

1. function numgrad = computeNumericalGradient(J, theta)  
2. % numgrad = computeNumericalGradient(J, theta)  
3. % theta: a vector of parameters  
4. % J: a function that outputs a real-number. Calling y = J(theta) will return the  
5. % function value at theta.   
6.     
7. % Initialize numgrad with zeros  
8. numgrad = zeros(size(theta));  
9.   
10. %% ---------- YOUR CODE HERE --------------------------------------  
11. % Instructions:   
12. % Implement numerical gradient checking, and return the result in numgrad.    
13. % (See Section 2.3 of the lecture notes.)  
14. % You should write code so that numgrad(i) is (the numerical approximation to) the   
15. % partial derivative of J with respect to the i-th input argument, evaluated at theta.    
16. % I.e., numgrad(i) should be the (approximately) the partial derivative of J with   
17. % respect to theta(i).  
18. %                  
19. % Hint: You will probably want to compute the elements of numgrad one at a time.   
20.   
21. epsilon = 1e-4;  
22.   
23. for i =1:length(numgrad)  
24.     oldT = theta(i);  
25.     theta(i)=oldT+epsilon;  
26.     pos = J(theta);  
27.     theta(i)=oldT-epsilon;  
28.     neg = J(theta);  
29.     numgrad(i) = (pos-neg)/(2*epsilon);  
30.     theta(i)=oldT;  
31.     if mod(i,100)==0  
32.        fprintf('Done with %d\n',i);  
33.     end;  
34. end;  
35.   
36.   
37. %% ---------------------------------------------------------------  
38. end  

cnnInitParams.m: 

[plain] view plain copy

1. function theta = cnnInitParams(imageDim,filterDim,numFilters,...  
2.                                 poolDim,numClasses)  
3. % Initialize parameters for a single layer convolutional neural  
4. % network followed by a softmax layer.  
5. %                              
6. % Parameters:  
7. %  imageDim   -  height/width of image  
8. %  filterDim  -  dimension of convolutional filter                              
9. %  numFilters -  number of convolutional filters  
10. %  poolDim    -  dimension of pooling area  
11. %  numClasses -  number of classes to predict  
12. %  
13. %  
14. % Returns:  
15. %  theta      -  unrolled parameter vector with initialized weights  
16.   
17. %% Initialize parameters randomly based on layer sizes.  
18. assert(filterDim < imageDim,'filterDim must be less that imageDim');  
19.   
20. Wc = 1e-1*randn(filterDim,filterDim,numFilters);  
21.   
22. outDim = imageDim - filterDim + 1; % dimension of convolved image  
23.   
24. % assume outDim is multiple of poolDim  
25. assert(mod(outDim,poolDim)==0,...  
26.        'poolDim must divide imageDim - filterDim + 1');  
27.   
28. outDim = outDim/poolDim;  
29. hiddenSize = outDim^2*numFilters;  
30.   
31. % we'll choose weights uniformly from the interval [-r, r]  
32. r  = sqrt(6) / sqrt(numClasses+hiddenSize+1);  
33. Wd = rand(numClasses, hiddenSize) * 2 * r - r;  
34.   
35. bc = zeros(numFilters, 1); %初始化为0  
36. bd = zeros(numClasses, 1);  
37.   
38. % Convert weights and bias gradients to the vector form.  
39. % This step will "unroll" (flatten and concatenate together) all   
40. % your parameters into a vector, which can then be used with minFunc.   
41. theta = [Wc(:) ; Wd(:) ; bc(:) ; bd(:)];  
42.   
43. end  

cnnParamsToStack.m: 

[plain] view plain copy

1. function [Wc, Wd, bc, bd] = cnnParamsToStack(theta,imageDim,filterDim,...  
2.                                  numFilters,poolDim,numClasses)  
3. % Converts unrolled parameters for a single layer convolutional neural  
4. % network followed by a softmax layer into structured weight  
5. % tensors/matrices and corresponding biases  
6. %                              
7. % Parameters:  
8. %  theta      -  unrolled parameter vectore  
9. %  imageDim   -  height/width of image  
10. %  filterDim  -  dimension of convolutional filter                              
11. %  numFilters -  number of convolutional filters  
12. %  poolDim    -  dimension of pooling area  
13. %  numClasses -  number of classes to predict  
14. %  
15. %  
16. % Returns:  
17. %  Wc      -  filterDim x filterDim x numFilters parameter matrix  
18. %  Wd      -  numClasses x hiddenSize parameter matrix, hiddenSize is  
19. %             calculated as numFilters*((imageDim-filterDim+1)/poolDim)^2   
20. %  bc      -  bias for convolution layer of size numFilters x 1  
21. %  bd      -  bias for dense layer of size hiddenSize x 1  
22.   
23.   
24. outDim = (imageDim - filterDim + 1)/poolDim;  
25. hiddenSize = outDim^2*numFilters;  
26.   
27.   
28. %% Reshape theta  
29. indS = 1;  
30. indE = filterDim^2*numFilters;  
31. Wc = reshape(theta(indS:indE),filterDim,filterDim,numFilters);  
32. indS = indE+1;  
33. indE = indE+hiddenSize*numClasses;  
34. Wd = reshape(theta(indS:indE),numClasses,hiddenSize);  
35. indS = indE+1;  
36. indE = indE+numFilters;  
37. bc = theta(indS:indE);  
38. bd = theta(indE+1:end);  
39.   
40.   
41. end  

　2013.12.30:

　　微博网友@路遥_机器学习：利用matlab自带的优化函数conv2,实现的mean-pooling，可以大大加快速度，大家可以参考。cnnPool.m文件里面：

[plain] view plain copy

1. tmp = conv2(convolvedFeatures(:,:,numFeature,numImage), ones(poolDim),'valid');   
2. pooledFeatures(:,:,numFeature,numImage) =1./(poolDim^2)*tmp(1:poolDim:end,1:poolDim:end);  

转自：http://www.cnblogs.com/tornadomeet/p/3468450.html