目的：打算使用深度學習的方式實現人臉關鍵點的檢測，第一步是要學習深度學習。

步驟：第一步在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