1: function patches = sampleIMAGES()

   2: % sampleIMAGES

   3: % Returns 10000 patches for training

   4: load IMAGES;    % load images from disk IMAGES是一个512*512*10的数组。

   5: patchsize = 8;  % we‘ll use 8x8 patches 

   6: numpatches = 10000;

   7: patches = zeros(patchsize*patchsize, numpatches);%这个patches是一个矩阵

   8: %每一列是一个64的列向量，是把每一个8*8的掩模变成这样一个向量，然后这个矩阵有

   9: %10000个列。也就是说从10幅512*512的图像中弄出来10000个掩模。现在初始化它。

  10: %% 从10幅512*512的图像中随机选出10000个8*8的掩模，然后制作 patches

  11: image_size=size(IMAGES);

  12: i=randi(image_size(1)-patchsize+1,1,numpatches);

  13: %在开区间（0，image_size(1)-patchsize+1）之间生成一个1行10000列的向量

  14: %也就是在一副512*512的图像中的行号中（不包括最后的7行，因为选了也构不成8*8的掩模了）

  15: %随意选出10000个数字组成向量A。下面这行是在列号中选10000个元素的向量B，

  16: %这样每一次从A和B中选出一个数a，b，然后a+7，b+7，这样a，b，a+7，b+7就组成了一个

  17: %8*8的掩模。（注意是选择哪个地方的掩模之前，要先从十幅图像中选一幅图像)

  18: %然后依次这么做，就得到了10000个掩模，然后把每个掩模弄成8*8的向量，

  19: %依次排到 patches的每一列上就得到了 patches。

  20: j=randi(image_size(2)-patchsize+1,1,numpatches);

  21: k=randi(image_size(3),1,numpatches);

  22: for num=1:numpatches

  23:     patches(:,num)=reshape(IMAGES(i(num):i(num)+patchsize-1,j(num):j(num)+patchsize-1,k(num)),1,patchsize*patchsize);

  24: end

  25: %% ---------------------------------------------------------------

  26: % For the autoencoder to work well we need to normalize the data

  27: % Specifically, since the output of the network is bounded between [0,1]

  28: % (due to the sigmoid activation function), we have to make sure 

  29: % the range of pixel values is also bounded between [0,1]

  30: patches = normalizeData(patches);

  31: display_network(patches(:,randi(size(patches,2),200,1)),8);%这个是随机从那10000个patch里面选出来200个patch，看看是不是

  32: %弄对了。

  33: end %这个end是说明patches这个function结束了，因为下面还有一个function，在一个文件里，所以一个function完了必须写end

  34:  

  35: %% ---------------------------------------------------------------

  36: function patches = normalizeData(patches)

  37: % Remove DC (mean of images). 

  38: patches = bsxfun(@minus, patches, mean(patches));

  39: % mean(patches)对patches每一列算出一个平均值。也就是每一组样本算出一个平均值。

  40: %mean(patches)是一个1*10000的行向量，这里 patches是一个64*10000的矩阵，按理说

  41: %不能相减，但这里 用到了bsxfun（@minus,），那么mean(patches)会自动复制64行，

  42: %然后也变成一个64*10000的矩阵，然后相减。

  43:  

  44: % Truncate to +/-3 standard deviations and scale to -1 to 1

  45: pstd = 3 * std(patches(:));%patches(:)是把这个矩阵变成一个向量。然后std，也就是

  46: %求所有10000个掩模所有像素值的标准差，这里如果你不是很明白为什么要把所以样本一块求标准差，

  47: %而不是一个样本组一个样本组的求标准差。并且还最后把所有像素的标准差*3。解释一下，以前求

  48: %数据都是一个样本组一个样本组进行归一化，那是把数据进行归一化，也就是均值为0，方差为1，但是

  49: %这样并不能保证所有的数据都在[-1,1]之间，只是方差为1。但是我们知道我们把所有数据减去均值

  50: %的话，那么所有数据的99.7%都会落在[-3*标准差，3*标准差]之间，所以我们只需要把剩下的0.03%的数据都

  51: %置成-3*标准差或3*标准差即可。这样所有数据都在[-3*标准差，3*标准差]之间，然后我们除以

  52: %3*标准差，那么素有数据都会在[-1,1]之间了。注意这里是将所有数据一块处理，而不需要一组一组样本处理，

  53: %因为这里并不需要处理后的每个样本满足标准正态分布:均值为0，方差为1那种关系，也不可能满足，所以

  54: %数据一块处理即可。

  55: patches = max(min(patches, pstd), -pstd) / pstd;

  56: %min(patches, pstd)意思是patches这个矩阵中元素如果大于pstd，那么这个元素就是prsd。

  57: %也就是上限是pstd，然后这里再max，也就是矩阵元素的下限是-petd。

  58:  

  59: % Rescale from [-1,1] to [0.1,0.9]

  60: patches = (patches + 1) * 0.4 + 0.1;

  61: end

   1: function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ...

   2:                                              lambda, sparsityParam, beta, data)

   3: %这个sparseAutoencoderCost函数的作用是：你给定神经网络的结构，网络参数值，L2正则项、稀疏约束项的权重，稀疏性参数。

   4: %那么这个函数就会给出这个神经网络的cost function的值，和cost function对各个参数的偏导。注意这个函数不优化各参数的

   5: %值，那个是优化minfun的活，这里是你给定这些参数，然后给出cost function 和偏导。如果你用优化算法比如梯度下降法优化参数

   6: %时，每一次迭代，你都要通过给出当下的参数，然后通过这个函数得到偏导，然后利用优化算法得到新的参数值，再带入这个函数

   7: %得到偏导，不停迭代，直到收敛。

   8: %注意visibleSize 应该是data的行数，既然data矩阵有了，就应该不需要再输入visibleSize了，确实

   9: %是这样，但是这里因为需要使用 visibleSize这个量，所以就又给出了一遍。你输入的visibleSize要保证等于

  10: %size（data，1）。

  11: % visibleSize: the number of input units (probably 64)

  12: % 输入层神经元个数，不包括始终为1的偏置节点

  13: % hiddenSize: the number of hidden units (probably 25) 这里神经元网络隐藏层是1层，

  14: %hiddenSize是这层神经元的个数,同样不包括始终为1的偏置节点

  15: % lambda: weight decay parameter L2正则项权重系数

  16: % sparsityParam: 就是稀疏性参数，一个标量，一般取为0.2等等

  17: % beta: weight of sparsity penalty term 稀疏项权重系数

  18: % data: Our 64x10000 matrix containing the training data.  So, data(:,i) is

  19: % the i-th training example. 样本数据  

  20: % The input theta is a vector (because minFunc expects the parameters to be a vector). 

  21: % We first convert theta to the (W1, W2, b1, b2) matrix/vector format, so that this 

  22: % follows the notation convention of the lecture notes. W1指第一层的参数W2第二层参数

  23: %输入的theta是个向量，向量元素的个数是神经元网络所有参数的个数这里面的元素怎么给参数，给哪些参数都行，

  24: %都可以按你的方式来，在这里是按这个顺序给参数的，theta的前

  25: %hiddenSize*visibleSize是神经网络第一层和W有关的参数，如果W下标的编号和

  26: %http://www.cnblogs.com/happylion/p/4193527.html第一个图中的w下标的编号一样，再结合下面reshape“按列”

  27: %进行重组成矩阵的该函数特性，这里把theta向量的前hiddenSize个元素给了W1矩阵的第一列，而W1矩阵的第一列各元素的含义是

  28: %输入层第一个神经元分别对应第二层隐藏层的各单元（不包括第二层偏置节点）的参数。W1矩阵的列数就是输入层神经元的个数。

  29: %这样theta向量里的元素的意义一旦都在神经元网络各参数中分配好了，之后训练计算等等之后，你就没法再把theta里的元素改

  30: %变网络参数的顺序或意义了。

  31: W1 = reshape(theta(1:hiddenSize*visibleSize), hiddenSize, visibleSize);%W1是一个hiddenSize*visibleSize的矩阵

  32: W2 = reshape(theta(hiddenSize*visibleSize+1:2*hiddenSize*visibleSize), visibleSize, hiddenSize);

  33: b1 = theta(2*hiddenSize*visibleSize+1:2*hiddenSize*visibleSize+hiddenSize);

  34: b2 = theta(2*hiddenSize*visibleSize+hiddenSize+1:end);

  35:  

  36: %这样神经网络的初始化就完成了。（注意这里初始化的代码表明只能有一层隐藏层）

  37: %% 初始化要求得的cost function 对各参数的偏导数

  38: % Cost and gradient variables (your code needs to compute these values). 

  39: % Here, we initialize them to zeros. 

  40: cost = 0;

  41: W1grad = zeros(size(W1)); %W1grad 也是一个矩阵，跟W1的行列数是一样的。里面的元素与神经网络里面的参数对应的情况和W1

  42: %里元素与神经网络参数对应的情况是一样的。第一列各元素的含义是cost function 对输入层第一个神经元分别对应第二层隐藏层的

  43: %各单元（不包括第二层偏置节点）的参数的偏导。

  44: W2grad = zeros(size(W2));

  45: b1grad = zeros(size(b1)); 

  46: b2grad = zeros(size(b2));

  47:  

  48: %% 求cost function 对各参数的偏导数

  49: %  Instructions: Compute the cost/optimization objective J_sparse(W,b) for the Sparse Autoencoder,

  50: %                and the corresponding gradients W1grad, W2grad, b1grad, b2grad.

  51: %

  52: % W1grad, W2grad, b1grad and b2grad should be computed using backpropagation.

  53: % Note that W1grad has the same dimensions as W1, b1grad has the same dimensions

  54: % as b1, etc.  Your code should set W1grad to be the partial derivative of J_sparse(W,b) with

  55: % respect to W1.  I.e., W1grad(i,j) should be the partial derivative of J_sparse(W,b) 

  56: % with respect to the input parameter W1(i,j).  Thus, W1grad should be equal to the term 

  57: % [(1/m) \Delta W^{(1)} + \lambda W^{(1)}] in the last block of pseudo-code in Section 2.2 

  58: % of the lecture notes (and similarly for W2grad, b1grad, b2grad).

  59: % 

  60: % Stated differently, if we were using batch gradient descent to optimize the parameters,

  61: % the gradient descent update to W1 would be W1 := W1 - alpha * W1grad, and similarly for W2, b1, b2. 

  62: % 

  63: %1.forward propagation

  64: data_size=size(data);

  65: biasPara_1=repmat(b1,1,data_size(2));%b1是第1层（输入层）总是1的那个偏置节点

  66: %与第二层各单元（不包括第二层偏置节点）对应的各个参数，b1这个向量元素的个数是第二层单元数-1个第二层偏置节点数

  67: %所以这里的biasPara_1这个矩阵是hiddenSize*样本个数的矩阵，并且矩阵每一列都是一样的。为什么要这么做，复制那么多列，

  68: %因为这里是所有样本一块处理的，不是一个样本一个样本的处理。因为data也是所有样本的数据，所以，你就能看出来下面的

  69: %W1*data+biasPara_1，所有样本是一块处理的：W1*data得到的矩阵第一列第一个元素是第一个样本各子元素在第二层

  70: %第一个元素上的线性叠加，所以还要加上b1的第1个元素；W1*data得到的矩阵第一列第2个元素是第一个样本各子元素在第2层

  71: %第1个元素上的线性叠加，所以还要加上b1的第2个元素；以此类推得到完整第一列，就是处理完第一个样本，然后处理第二列，

  72: %就是处理第二个样本，以此类推..所以W1*data+biasPara_1这个矩阵第m列第n行的元素就是第m个样本，对第二层第n单元

  73: %（不包括偏置节点，因为偏置节点不需要什么sigmoid函数，也不需要有输入，它始终是1，可以认为激活值始终是1）

  74: %的sigmoid函数的输入Z。

  75: biasPara_2=repmat(b2,1,data_size(2));

  76: active_value2=sigmoid(W1*data+biasPara_1);

  77: active_value3=sigmoid(W2*active_value2+biasPara_2);

  78: %W2*active_value2+biasPara_2这个矩阵第m列第n行的元素就是第m个样本，对第3层第n单元的sigmoid函数的输入Z。

  79: %active_value2这个矩阵每一列是每一个样本前向算法中第二层各单元的激活值（不包括偏置节点，因为偏置节点

  80: %相当于激活值始终为1）。

  81: %active_value3是第三层各单元的激活值。

  82:  

  83: %2.computing cost function 

  84: ave_square=sum(sum((active_value3-data).^2)./2)/data_size(2);

  85: %这是cost

  86: %function的均方误差项，因为这里data的label还相当于是data（自编码），首先说一下一个矩阵“.^2”是矩阵里每个元素都平方。

  87: %sum((active_value3-data).^2）得到的是一个行向量，每一个元素是每一个样本通过神经网络得到的输出层单元激活值（是一个向量）

  88: %与每一个样本的label（也就是它自身）（也是一个向量）的误差平方（计算是通过两个向量对应的各个元素分别对应做差然后平方和）

  89: %然后sum((active_value3-data).^2）得到的这个向量各元素都乘以1/2（也就是每个样本的误差平方乘以一个系数1/2），然后再sum

  90: %一下，就是把所有样本的误差平方累加起来，最后再除以样本的个数data_size(2)，就是cost function的均方误差项。

  91: weight_decay=lambda/2*(sum(sum(W1.^2))+sum(sum(W2.^2)));

  92: %这是const function 中的L2正则项

  93: active_average=sum(active_value2,2)./data_size(2);%sum(active_value2,2)这里是active_value2按行相加，也就是对隐藏层每一个

  94: %神经元，将所有样本关于这个神经元的激活值累加起来。sum（，2）是按行相加，如果只是sum（），就是默认的按列相加。matlab中

  95: %很多函数默认的都是按列进行某种操作。因此这里得到的是一个列向量。然后点除data_size(2)就是向量每个元素都除以样本个数

  96: %这样得到的列向量active_average的每个元素就是隐藏层每一个神经元对于所有样本的平均激活度。

  97: p_para=repmat(sparsityParam,hiddenSize,1);

  98: %因为下面那行代码中有p_para./active_average所以要弄成和active_average一样size的列向量，才能与active_average里的元素点除。

  99: sparsity=beta.*sum(p_para.*log(p_para./active_average)+(1-p_para).*log((1-p_para)./(1-active_average)));

 100: %这是稀疏项

 101: cost=ave_square+weight_decay+sparsity;%cost function就完成了。

 102:  

 103: %求残差

 104: delta3=(active_value3-data).*(active_value3).*(1-active_value3);

 105: %首先说博客中求残差的公式是按一个样本的一个神经元来定位的，而这里用matlab求的时候是所有的样本（一列是一个样本），

 106: %一层所有的神经元一块运算的，这是matlab 矩阵运算的特点，避免了for循环。所以这里跟公式中略有不同。虽然是一起运算的，

 107: %但每一个神经元残差的计算还是要根据公式来，所以这里用到了点乘：矩阵对应元素之间运算。

 108: %公式在http://www.cnblogs.com/happylion/p/4193527.html可以看到。

 109: %这里得到的delta3是一个矩阵，每一列是一个样本输入神经网络后，第三层各个神经元的残差，列数代表样本数。

 110: active_average_repmat=repmat(sum(active_value2,2)./data_size(2),1,data_size(2));

 111: %sum(active_value2,2)./data_size(2)就是active_average，然后repmat，因为active_average是一个列向量，再平铺

 112: %data_size(2)列，那么active_average_repmat就是一个矩阵，弄成矩阵的的形式是为了后面所有样本第二层所有神经元

 113: %一起运算。

 114: default_sparsity=repmat(sparsityParam,hiddenSize,data_size(2));

 115: %default_sparsity也是一个hiddenSize*data_size(2)矩阵，弄成这个形式也是为了下面所有样本第二层所有神经元一起的矩阵

 116: %运算。

 117: sparsity_penalty=beta.*(-(default_sparsity./active_average_repmat)+((1-default_sparsity)./(1-active_average_repmat)));

 118: delta2=(W2‘*delta3+sparsity_penalty).*((active_value2).*(1-active_value2));

 119: %求第二层的残差的时候跟输出层略有不同，因为第二层是隐藏层，又因为稀疏项中是关于隐藏层的内容，所以这里把稀疏项的内容

 120: %融入到残差里面。这样由残差计算偏导的公式就不用改变了。同样这里是所有样本第二层所有神经元的一起计算的矩阵计算，

 121: %所以也用到了点乘，矩阵之间对应元素的运算就是博客中的公式。((active_value2).*(1-active_value2))这一项是f‘(z(i))

 122: %W2‘*delta3得到的是一个矩阵，矩阵的每一列是一个样本输入神经网络后，第二层各个神经元的残差公式的第一项部分（可以看

 123: %http://www.cnblogs.com/happylion/p/4193527.html中第二个图就会很明白了），列数代表样本数。

 124:  

 125: %求各参数偏导

 126: W2grad=delta3*active_value2‘./data_size(2)+lambda.*W2;

 127: %delta3*active_value2‘得到的是一个矩阵，矩阵的第一个元素是所有样本对应的 小cost function

 128: %对第二层的W11的偏导 的累加形成的总cost function对第二层的W11的偏导。矩阵第一行第二个元素是

 129: %总cost function对第二层的W21的偏导。所以你看这里的W下标的顺序和前面的W矩阵是一样的。然后这个矩阵的元素都除以样本

 130: %个数，得到的就是cost function中只有均方误差项时，对第二层各个参数的偏导。然后这个矩阵每个元素加上lambda倍的W2这个

 131: %矩阵对应的元素，得到的W2grad矩阵的元素是总cost function对第二层各参数的偏导。

 132: W1grad=delta2*data‘./data_size(2)+lambda.*W1;

 133: b2grad=sum(delta3,2)./data_size(2);%b2grad向量是总cost function对第二层偏置节点各参数的偏导

 134: b1grad=sum(delta2,2)./data_size(2);

 135:  

 136: %-------------------------------------------------------------------

 137: % After computing the cost and gradient, we will convert the gradients back

 138: % to a vector format (suitable for minFunc).  Specifically, we will unroll

 139: % your gradient matrices into a vector.

 140:  

 141: grad = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)];%正好矩阵（:）这种形式变成一个向量也是按列的，正好和

 142: %reshape（）按列的方式一样。所以元素与网络参数对应的顺序没有变化。

 143:  

 144: end

 145:  

 146: %-------------------------------------------------------------------

 147: % Here‘s an implementation of the sigmoid function, which you may find useful

 148: % in your computation of the costs and the gradients.  This inputs a (row or

 149: % column) vector (say (z1, z2, z3)) and returns (f(z1), f(z2), f(z3)). 

 150: %若果x是个矩阵，也是可以的，igmoid(x)也是一个矩阵，里面的每个元素就是对应x矩阵对应的元素进行sigmoid计算。

 152: function sigm = sigmoid(x)

 154:     sigm = 1 ./ (1 + exp(-x));

 155: end

   1: function numgrad = computeNumericalGradient(J, theta)

   2: %% computeNumericalGradient函数的功能是：给定一个方程J(theta)，给定的方式是@你想使用的方程，

   3: %然后你再给一个theta的值，那么这个函数就可以利用求导定义的那个公式来求出这个J方程在theta这点的导数，如果

   4: %theta是一个向量，那么就是求出J方程在theta这点所有未知量的偏导数，因此得出的结果是一个向量，以 numgrad

   5: %给出。

   6: numgrad = zeros(size(theta));%初始化

   7: EPSILON=0.0001;

   8: for i=1:size(theta)  %每一个i是J方程对theta的第i个变量求偏导数。

   9:     theta_plus=theta;

  10:     theta_minu=theta;

  11:     theta_plus(i)=theta_plus(i)+EPSILON;%用定义公式求偏导，那么theta这点的很小的增量或减量中只有第i个

  12:     %变量增加或减少了EPSILON

  13:     theta_minu(i)=theta_minu(i)-EPSILON;

  14:     numgrad(i)=(J(theta_plus)-J(theta_minu))/(2*EPSILON);

  15: end

  16: end

   1: %% 当我们检查sparseAutoencoderCost函数写的对不对或这个函数得到的偏导数是不是跟导数的定义一致，

   2: %我们随机给定一个theta 即可，如果这个随机给出的theta使上述两个函数的偏导数求得的一致，那么我们就认为任何theta

   3: %都会使之相应的一致。并且我们不必把所有样本都拿来检验，只需要把其中的10个做为样本，把隐藏层单元弄成2，3个即可。

   4: %但输入层单元个数不能改变，因为它等于样本矩阵的行数。最后我们用norm(numgrad-grad)/norm(numgrad+grad)

   5: %进行判定numgrad、grad两个向量的相似程度，如果结果 <1e-9,那么就说他们相似到我们希望的程度了，即检查结果满意！

   6: visibleSize = 8*8;   

   7: hiddenSize = 3;    

   8: sparsityParam = 0.01;   

   9: lambda = 0.0001;   

  10: beta = 3;          

  11: patches = sampleIMAGES;

  12: theta = initializeParameters(hiddenSize, visibleSize);

  13: [cost, grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, lambda, ...

  14:                                      sparsityParam, beta, patches(:,1:10));

  15: numgrad = computeNumericalGradient( @(x) sparseAutoencoderCost(x, visibleSize, ...

  16:                                                   hiddenSize, lambda, ...

  17:                                                   sparsityParam, beta, ...

  18:                                                   patches(:,1:10)), theta);

  19: disp([numgrad grad]); 

  20: diff = norm(numgrad-grad)/norm(numgrad+grad);

  21: disp(diff); 

   1: visibleSize = 8*8;   

   2: hiddenSize = 25;    

   3: sparsityParam = 0.01;   

   4: lambda = 0.0001;   

   5: beta = 3;          

   6: patches = sampleIMAGES;

   7: theta = initializeParameters(hiddenSize, visibleSize);

   8: %  Use minFunc to minimize the function

   9: addpath minFunc/

  10: options.Method = ‘lbfgs‘; % Here, we use L-BFGS to optimize our cost

  11:                           % function. Generally, for minFunc to work, you

  12:                           % need a function pointer with two outputs: the

  13:                           % function value and the gradient. In our problem,

  14:                           % sparseAutoencoderCost.m satisfies this.

  15: options.maxIter = 400;      % Maximum number of iterations of L-BFGS to run 

  16: options.display = ‘on‘;

  17: [opttheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ...

  18:                                    visibleSize, hiddenSize, ...

  19:                                    lambda, sparsityParam, ...

  20:                                    beta, patches), ...

  21:                               theta, options);

  22: %因为每次迭代在sparseAutoencoderCost函数只需要更改参数网络，网络结构不变，所以用到了

  23: %@(p) sparseAutoencoderCost(p, visibleSize, hiddenSize, lambda, sparsityParam, beta, patches), theta, options)