ufldl 深度学习入门 第一发：基于BP网络实现稀疏自编码器

目的：打算使用深度学习的方式实现人脸关键点的检测，第一步是要学习深度学习。

步骤：第一步在ufldl上面学习深度学习的算法基础知识，然后找博客上基于python调用theano库实现人脸关键点检测的算法，看懂后基于C++实现，然后用java实现app，调用C++实现的算法。

ufldl上的第一章是基于BP网络实现稀疏自编码器，在matlab上实现。

稀疏自编码器的实现：ufldl上已经给出了整体的框架，我们需要将三个地方补全。由于我对matlab各种函数不熟悉，所以采用的是参考别人实现的代码，然后去理解的方式，等到ufldl教程后面几课越来越熟悉，再自己来实现。

首先是补全sampleIMAGES.m中的代码，该代码从IMAGES.mat中随机提取8×8×10000张sample patches，用来做输入样本。

    tic  
    image_size=size(IMAGES);  
    i=randi(image_size(1)-patchsize+1,1,numpatches);   %产生1*10000个随机数 范围在[1,image_size(1)-patchsize+1]之间  
    j=randi(image_size(2)-patchsize+1,1,numpatches);  
    k=randi(image_size(3),1,numpatches);              % 随机的选取图片 10000次  
    for num=1:numpatches  
            patches(:,num)=reshape(IMAGES(i(num):i(num)+patchsize-1,j(num):j(num)+patchsize-1,k(num)),1,patchsize*patchsize);  
    end  
    toc  

说明：

tic和toc用来计时，返回tic-toc之间代码运行花费的时间;

IMAGES是512×512×10的数组，所以size（IMAGES）返回512 512 10的向量;

i=randi(512-8+1,1,10000)，会返回一个1×10000的随机数组，数组中的元素属于（0，505]，也就是[1，505];（ps：上面介绍randi函数说产生的是开区间，但是我试了发现并不是开区间，是可以取到505的，不过这里并没有什么影响，就不去纠结这一点了）

所以i j k都是1×10000的数组，元素大小位于1-505;

sample=IMAGES(i(num):i(num)+patchsize-1,j(num):j(num)+patchsize-1,k(num))

上面这句话采用全下标的方式访问3维数组IMAGES的元素，行从i（num）到i（num）+7，列从j（num）到j（num）+7，页选择k（num）;

reshape（sample，1，64）将sample变成了1×64的行向量，然后赋值给了patches;

sample到图片后又将图片的数据全部归一化到了[0.1，0.9]：patches = normalizeData(patches)

这里的归一化函数，之后再分析，之所以要归一化到[0.1，0.9]，是因为激活函数sigmoid函数的输出范围是[0，1]，所以需要将输入归一化到函数的输出范围内，这样才能训练。

其次是实现sparseAutoencoder.m，稀疏自编码

%1.forward propagation
data_size=size(data);
              %获得data的维数信息，64×10000
active_value2=repmat(b1,1,data_size(2));     
               %扩展b1，得到结果为25×10000，之前b1为25×1，对应的是一个patch，现在有10000个patch，相应的变成10000columns
active_value3=repmat(b2,1,data_size(2));    
               %同上，对应b2
active_value2=sigmoid(W1*data+active_value2);     
               %计算中间隐藏层的激活值，z=w1×data+b1 对应25×10000
active_value3=sigmoid(W2*active_value2+active_value3);     
               %计算输出层的激活值，z=w2×a2+b2 对应64×10000
%2.computing error term and cost
ave_square=sum(sum((active_value3-data).^2)./2)/data_size(2);     
               %计算均方误差
weight_decay=lambda/2*(sum(sum(W1.^2))+sum(sum(W2.^2)));    
               %计算权重衰减项

p_real=sum(active_value2,2)./data_size(2);     
               %计算平均活跃度，sum（x，2）将每行相加，结果为25×1的列向量
p_para=repmat(sparsityParam,hiddenSize,1);    
               %将参数p，repmat为25×1的列向量
sparsity=beta.*sum(p_para.*log(p_para./p_real)+(1-p_para).*log((1-p_para)./(1-p_real)));     %求得稀疏性限制项
cost=ave_square+weight_decay+sparsity;     
               %3项求和得到总的代价函数，但是为什么要求总的代价函数呢？
               %貌似并不需要啊，只需要求偏导就可以了啊，确实这4句话并非必须
               %其实是因为后面的computeNumercialGradient函数，需要总的代价cost
               %这样就可以通过导数的定义来计算gradient了
delta3=(active_value3-data).*(active_value3).*(1-active_value3);     
               %计算输出层的残差，结果形式是64×10000
average_sparsity=repmat(sum(active_value2,2)./data_size(2),1,data_size(2));     
               %计算平均活跃度，并且repmat成25×10000的矩阵形式
default_sparsity=repmat(sparsityParam,hiddenSize,data_size(2));     
               %把默认参数p，repmat成计算需要的25×10000矩阵形式
sparsity_penalty=beta.*(-(default_sparsity./average_sparsity)+((1-default_sparsity)./(1-average_sparsity)));  
               %计算 计算delta时需要的稀疏惩罚项
delta2=(W2'*delta3+sparsity_penalty).*((active_value2).*(1-active_value2));     
               %计算中间隐藏层的残差，这里加入了稀疏惩罚项
%3.backword propagation 后向传播，更新w和b的值
W2grad=delta3*active_value2'./data_size(2)+lambda.*W2;       
               %64×10000  ×   10000×25 然后求平均值得到64×25的W2grad
W1grad=delta2*data'./data_size(2)+lambda.*W1;        
               %25×10000  ×  10000×64 然后求平均值得到25×64的W1grad
b2grad=sum(delta3,2)./data_size(2);      
               %得到64×1
b1grad=sum(delta2,2)./data_size(2);     
               %得到25×1
               %可是问题是这里的代码并没有对w b进行更新迭代啊，
               %只是求了一遍，更新迭代的代码在哪里？ 
               %andrew ng 给出了后面要用的迭代求解的代码，会重复调用这个函数
               % minFunc函数就是通过LBFGS来快速迭代求解的 

最后数值计算导数，用作梯度检验，computeNumercialGradient.m

EPSILON=0.0001;
for i=1:size(theta)
               % for循环，从1到3289（25×64+64×25+25+64）
    theta_plus=theta;
    theta_minu=theta;
               % 将theta赋值给两个计算中用到的变量
    theta_plus(i)=theta_plus(i)+EPSILON;
    theta_minu(i)=theta_minu(i)-EPSILON;
                % 将原有w1 w2 b1 b2 中的某个值做一个很小的更改
                % 通过更改后的cost的变化，根据导数的定义计算
    numgrad(i)=(J(theta_plus)-J(theta_minu))/(2*EPSILON);
                % j是一个函数句柄变量，可以用来调用函数sparseAutoencoderCost
                % numgrad is 3289×1 vector
                % 为什么numgrad是3289×1的vector，J()的输出是[cost，grad]，这里还没有深究？？
end

基于上面的三段程序，加上andrew ng提供的程序，就可以完成稀疏自编码器的；

注意梯度检验这一段程序，用来检验前面写的程序是否正确，在确认正确后，train时，将该段程序注释掉

这样才能比较快的完成train，否则梯度检验这一段程序很慢的；

最后给出训练后的结果：

矢量化编程：

由于前面代码的实现已经是矢量化了，所以不需要大的更改，只需要改变读取输入数据的方式即可。

将STEP1 Implement sampleIMAGES 改成如下即可

images=loadMNISTImages('train-images-idx3-ubyte');   % image is 784*60000 matrix, 784=28*28
display_network(images(:,1:100));                    % Show the first 100 images
patches = images(:,1:10000);

最后得到的训练结果如图：

matlab中涉及到的函数如下，并附上使用简介：

1 cumtrapz(a);计算数组a的数值积分，比如a=[1 2 3 4 ]，cumtrapz（a）=[0 1.5 4 7.5 ]，因为a对应的图形的面积在这4个点处分别是0  1.5  4  7.5，默认a的每个元素之间的间距是1。

2 ndims（a）：返回数组a的维数=2。a=[1 2 3 4 ]或者a=1，a的维数都是2，有行和列。

3 size（a）：返回数组a的所有维度的值，所以返回的是一个向量值。比如a=a=[1 2 3 4 ]，a有两个维度，第一维度是行维度=1，第二维度是列维度=4。

4 随机函数

rand 生成均匀分布的伪随机数。分布在（0~1）之间

 主要语法：rand(m,n)生成m行n列的均匀分布的伪随机数
          rand(m,n,'double')生成指定精度的均匀分布的伪随机数，参数还可以是'single'
          rand(RandStream,m,n)利用指定的RandStream(我理解为随机种子)生成伪随机数
randn 生成标准正态分布的伪随机数（均值为0，方差为1）
   主要语法：和上面一样
randi 生成均匀分布的伪随机整数
  主要语法：randi（iMax）在开区间（0，iMax）生成均匀分布的伪随机整数
          randi（iMax，m，n）在开区间（0，iMax）生成mXn型随机矩阵
           r =randi([iMin,iMax],m,n)在开区间（iMin，iMax）生成mXn型随机矩阵

5 reshape（a，2，6），a=[1 2 3;4 5 6;7 8 9;10 11 12]，reshape（a，2，6）=[1 7 2 8 3 9；4 10 5 11 6 12]，具体的原理百度；按列拼接，然后抽取元素，组成2组行向量。

6 sum函数，sum(x,2)表示矩阵x的横向相加，求每行的和，结果是列向量；而缺省的sum(x)就是竖向相加，求每列的和，结果是行向量
7 disp（）函数，用于显示数组。a=‘hello’，disp（a）=hello；a=[1 2]，b=[3 4]，disp（[a，b]）=1 2 3 4