【deep learning精華部分】稀疏自編碼提取高階特征、多層微調完全解釋及代碼逐行詳解

我們前面已經講了如何訓練稀疏自編碼神經網絡，當我們訓練好這個神經網絡后，當有新的樣本輸入到這個訓練好的稀疏自編碼器中后，那么隱藏層各單元的激活值組成的向量就可以代表（因為根據稀疏自編碼，我們可以用來恢復），也就是說就是在新的特征下的特征值。每一個特征是使某一個取最大值的輸入。假設隱藏層單元有200個，那么就一共有200個特征，所以新的特征向量有200維。特征顯示情況在前面博客中已經給出，我們把這時候的特征稱為一階特征。

我們知道腦神經在處理問題，比如看一個圖片的時候，也不只使用了一層的神經，而是使用了多層的神經元不斷的將看到的圖片內容抽象，（更詳細的可以看這里） 。因為每一層都是一層非線性的疊加，所以每增加一層隱藏層，那么神經網絡的表達能力就增強一步，所以我們將要做的是提取二階特征。

注意我們這里要得到二階特征並不是使用像如下圖所示的“增加版”稀疏自編碼神經網絡（在原基礎增加一個隱藏層）：

我們為什么不適用這樣的一個網絡呢？原因是這樣的（引自這里）：

‘’

雖然幾十年前人們就發現了深度網絡在理論上的簡潔性和較強的表達能力，但是直到最近，研究者們也沒有在訓練深度網絡方面取得多少進步。 問題原因在於研究者們主要使用的學習算法是：首先隨機初始化深度網絡的權重，然后使用有監督的目標函數在有標簽的訓練集 上進行訓練。例如通過使用梯度下降法來降低訓練誤差。然而，這種方法通常不是十分奏效。這其中有如下幾方面原因：

數據獲取問題

使用上面提到的方法，我們需要依賴於有標簽的數據才能進行訓練。然而有標簽的數據通常是稀缺的，因此對於許多問題，我們很難獲得足夠多的樣本來擬合一個復雜模型的參數。例如，考慮到深度網絡具有強大的表達能力，在不充足的數據上進行訓練將會導致過擬合。

局部極值問題

使用監督學習方法來對淺層網絡（只有一個隱藏層）進行訓練通常能夠使參數收斂到合理的范圍內。但是當用這種方法來訓練深度網絡的時候，並不能取得很好的效果。特別的，使用監督學習方法訓練神經網絡時，通常會涉及到求解一個高度非凸的優化問題（例如最小化訓練誤差 ，其中參數 是要優化的參數。對深度網絡而言，這種非凸優化問題的搜索區域中充斥着大量“壞”的局部極值，因而使用梯度下降法（或者像共軛梯度下降法，L-BFGS等方法）效果並不好。

梯度彌散問題

梯度下降法（以及相關的L-BFGS算法等）在使用隨機初始化權重的深度網絡上效果不好的技術原因是：梯度會變得非常小。具體而言，當使用反向傳播方法計算導數的時候，隨着網絡的深度的增加，反向傳播的梯度（從輸出層到網絡的最初幾層）的幅度值會急劇地減小。結果就造成了整體的損失函數相對於最初幾層的權重的導數非常小。這樣，當使用梯度下降法的時候，最初幾層的權重變化非常緩慢，以至於它們不能夠從樣本中進行有效的學習。這種問題通常被稱為“梯度的彌散”.

與梯度彌散問題緊密相關的問題是：當神經網絡中的最后幾層含有足夠數量神經元的時候，可能單獨這幾層就足以對有標簽數據進行建模，而不用最初幾層的幫助。因此，對所有層都使用隨機初始化的方法訓練得到的整個網絡的性能將會與訓練得到的淺層網絡（僅由深度網絡的最后幾層組成的淺層網絡）的性能相似。

                                                                                                                                                                        ‘’

上面這個解釋不光單指稀疏自編碼神經網絡（輸出等於）的各種困難，所有的如上圖那樣直接增加隱藏層，然后利用

樣本訓練都會出現上面的困難。那么應該怎么辦呢？

我們這里使用了一種逐層貪婪訓練的策略。簡單來說，逐層貪婪算法的主要思路是每次只訓練網絡中的一層，即我們首先訓練一個只含一個隱藏層的網絡，跟前面方法一樣，我們根據輸出層與樣本標記的誤差平方作為cost function，然后使用梯度下降法或L-BFGS等方法（求梯度時結合前向后向傳播算法）最小化cost function得到第一隱藏層這一層的參數（第一隱藏層這一層指的是和的參數（書寫方便，這里和下面都只用代替，但是記得是包括的），注意訓練好以后我們只需要如下圖的參數，的參數在之后訓練新增加的隱藏層或最后的微調中是不需要的）。如圖所示（樣本的label是一個二維的向量，隱藏層用了三個神經元）：

如果是稀疏自編碼，其他都一樣，只是把前面提到的label，即樣本的標記換成了樣本本身，那么網絡輸出的就是原樣本，如圖所示（隱藏層用了四個神經元）：

 

然后當這層網絡訓練結束之后才開始訓練第二個隱藏層的網絡，訓練是這樣的：訓練的第二個隱藏層神經網絡的輸入是第一個隱藏層神經網絡隱藏層的激活值，輸出還是原樣本的label，然后跟上面一樣訓練第二個隱藏層的網絡，同樣地，對於這個神經網絡，我們需要的也是的參數（包括b（1））（為了與上面的區分，這里寫成），的參數在之后訓練新增加的隱藏層或最后的微調中是不需要的,如圖所示（隱藏層用了三個神經元）：

 

如果是稀疏自編碼，其他都一樣，只是把前面提到的樣本標記label換成了這一層神經網絡的輸入：，你可能發現不同了，如果是監督學習的時候，我們進行每一層貪婪訓練時，每一個網絡的輸出都是一樣的，都是原樣本的label。而如果每一層都是稀疏自編碼時，我們每一個網絡的輸出是相互不一樣的，都等於本層訓練網絡的輸入。如圖所示（隱藏層用了三個神經元）：

之后如果還需要增加層次，跟前面一樣，只是如果是稀疏自編碼的話，那么最終得到的是在新特征下的特征值，所以跟監督學習不一樣，整個過程並沒有分類，所以我們還需要增加一個分類器，比如softmax回歸，來對非監督學習到的新特征值進行分類，具體地：

當我們用一些不含label的樣本完成了稀疏自編碼逐層貪婪訓練后，也就是我們找到了如果用新的特征去表示原數據，這時候我們把另外一些含label的樣本，通過下圖所示的神經網絡，那么得到的就是這些樣本在之前學習到的特征下的特征值。注意這個網絡的參數就是前面逐層訓練時得到的和。那么在之后softmax訓練時，我們的訓練樣本和測試樣本就是這些新的特征值（即用這些新得特征值來表示樣本）：

……………………………………………………………………………………………………………………..圖1

得到新特征下的訓練樣本和測試樣本后，我們接下來要做的跟以前softmax訓練就一樣了，訓練得到softmax參數θ。

 

可能讀者覺得這樣基本就可以說ok了，其實這還遠遠不夠，我們現在得到的參數，只是逐層貪婪訓練得到的稀疏自編碼參數和最后一層監督訓練的softmax參數，但是我們真正想要的是那種純監督學習對整個神經網絡一塊訓練后的參數（雖然根據前面提到的缺陷，我們很難得到，但是我們還是想得到的這樣的參數），所以我們剛才做的工作其實相當於把整個神經網絡的參數由隨機的參數初始化變成了由‘逐層貪婪訓練得到稀疏自編碼參數和監督訓練softmax回歸得到最后一層參數’這種方法的參數初始化，我們這樣做的目的就是期待這樣的參數初始化方法可以（比參數隨機初始化）使 cost function最小化時更容易收斂。所以我們的任務還沒有結束，並且接下來的過程能使網絡性能大大提高，這一步我們俗稱微調（fine-tuning）。

上述如果想做微調，首先需要注意的一點是，我們訓練softmax回歸模型的樣本要和逐層訓練稀疏編碼的樣本都是來着同一組樣本，也就是說圖1輸入的還是原來訓練稀疏編碼的那些樣本（如果你用新的樣本的話，就代表之前訓練的參數已經訓練好了，所以還是要用原來的樣本，這才能進行微調），最后得到的去訓練和測試softmax回歸模型得到圖2最后的輸出p(y=1|x)、p(y=2|x)、p(y=3|x)，即我們把這整個過程看成一個整體，如下圖2所示；我們利用一組樣本訓練出上面所說的W(1)、WW(1)和θ后，然后我們把原樣本一個個通過下面這個神經網絡中，在n-1層的神經元的激活值是各個樣本對應的，然后通過最后一層softmax模型得到最后的輸出。然后根據樣本的label y（i）和這個神經網絡的輸出p(y=1|x)、p(y=2|x)、p(y=3|x)，形成cost function，整個模型的cost function為：，

 

在對cost function求最小值的過程中，我們將cost function對各個參數求偏導，然后利用如梯度下降法迭代得到每一次迭代后的參數，最后迭代到cost function收斂后最終參數。在這里我們要注意一個地方：我們每一次迭代中求cost function 對前n-1層的網絡參數的導數是根據前向后向傳播算法計算的，而cost function 對最后的softmax層參數的導數我們是對上面的J（θ）直接求的，像以前博客中求softmax模型參數導數一樣。對於求前面幾層參數的導數，因為前面n-1行的網絡就是前面《神經網絡中的參數的求解：前向和反向傳播算法》博客中的網絡，前面后向傳播算法等大體也是上面的步驟。但是因為有最后softmax層，我們模型的cost function形式不一樣，所以還是多少有些出入，這里我們把不同的地方重新求解一下：

我們進行后向算法的時候，只對圖2中前n-1層進行后向傳播，求前n-1層神經元的殘差。

現在我們先對圖2中神經元求殘差，從上面提到的博客中我們知道是，f為sigmoid函數。

我們在利用前向后向算法時，是一個樣本一個樣本進行的，所以每一個樣本下模型的’小cost function‘是，i指的是第i個樣本下的各個值。下面我們’小cost function’對求導，因為是一樣向量，對它求導就是對向量里面的所以元素分別求導，我們現在由形成的’小 cost function’對的第m個元素求導：

 

這里的θ1m指的是前面博客定義的θ1（博客中定義的θ矩陣每一行分別為θ1轉置，θ2轉置，，，θk轉置）向量的第m個元素。

然后小 cost function’對所有的元素求導的向量形式就可以表示為：。I,P都是當前下的列向量（也就是在一個樣本下）。這是一個列向量，里面的每一個元素就是小 cost function’對的某個元素求導的值。以上這些是某一個下的計算，我們要對所有的分別進行如上操作（在matlab中采用矩陣（矢量）運算中I,P是一個矩陣，每一列是對應的樣本下的情況，那么也是一個矩陣，每一列對應一個樣本的情況）。

之后的步驟跟《神經網絡中的參數的求解：前向和反向傳播算法》博客中步驟是一樣的了，這里就不再重復了（還是有一點不同之處，這里我們說的cost function都是在沒有稀疏項和L2正則項的，我們要加入這些正則項。而在以前的稀疏編碼中只有一層隱藏層，而這里是有兩層，所以在代碼中還是有些不同的，具體的看代碼注釋）。

 

 

………………………………………………………………………………………………………………….圖2

 

 

上述這種方法稱為逐層貪婪訓練的半監督的機器學習方法（又叫自學習方法）：原始樣本提取新特征階段利用稀疏自編碼進行無監督學習，但學習到一定程度的特征后，然后利用監督學習方法（如softmax回歸）進行訓練分類等。如果是上面提到的第一類直接進行有監督的學習，即每層進行貪婪訓練時使用樣本的label作為每一層訓練網絡的輸出，那么就不需要最后的softmax等的分類學習了，因為它一開始就進行了監督學習，（這可能也是上面黃色字體所表達的不同的原因吧）所以它稱為逐層貪婪訓練的監督學習。現在比較常用的是前者：逐層貪婪訓練的半監督學習。

那么這里的逐層貪婪訓練的半監督學習跟最前面提到的“增加版”稀疏自編碼神經網絡（也就是在單隱藏層增加n個隱藏層后整體直接進行監督學習，也叫純監督學習、朴素監督學習或直接叫監督學習）相比，有哪些優點呢（引自這里）：

‘’

數據獲取

雖然獲取有標簽數據的代價是昂貴的，但獲取大量的無標簽數據是容易的。自學習方法（self-taught learning）的潛力在於它能通過使用大量的無標簽數據來學習到更好的模型。具體而言，該方法使用無標簽數據來學習得到所有層（不包括用於預測標簽的最終分類層） 的最佳初始權重。相比純監督學習方法，這種自學習方法能夠利用多得多的數據，並且能夠學習和發現數據中存在的模式。因此該方法通常能夠提高分類器的性能。

更好的局部極值

當用無標簽數據訓練完網絡后，相比於隨機初始化而言，各層初始權重會位於參數空間中較好的位置上。然后我們可以從這些位置出發進一步微調權重。從經驗上來說，以這些位置為起點開始梯度下降更有可能收斂到比較好的局部極值點，這是因為無標簽數據已經提供了大量輸入數據中包含的模式的先驗信息。

現在看一下程序：

我們還是去識別字體，數據庫都和以前的一樣。

首先我們先給出微調的這整個網絡參數：

%%  Here we provide the relevant parameters values that will
%  allow your sparse autoencoder to get good filters;
 clear all ; clc
inputSize = 28 * 28;%輸入層的神經元個數（不包括+1的截距神經元，下面的hiddenSizeL1，2也不包括，因為在進行稀疏自編碼的時候，以前我們寫的函數
%使我們只需要輸入不包括+1截距神經元信息的神經元個數和樣本即可。但是下面的softmax層，我們在輸入的時候我們就要添加+1層截距神經元信息，這點要
%注意，不然很容易出bug）
numClasses = 10;       %類別有10個，label分別為1，2....10
hiddenSizeL1 = 200;    % Layer 1 Hidden Size  如果你想展示你的二階特征，這個值必須是一個數平方的值
hiddenSizeL2 = 200;    % Layer 2 Hidden Size
sparsityParam = 0.1;   % desired average activation of the hidden units.
                       % (This was denoted by the Greek alphabet rho, which looks like a lower-case "p",
                       %  in the lecture notes).  兩層隱藏層的各神經元的平均激活度都是這個值
lambda = 3e-3;         % 稀疏自編碼各層的L2范數權重  這個 3e-3參數值如果效果不好，可以自己再調一調
lambda2 = 1e-4;        %softmax層L2范數權重 這個1e-4 參數值如果效果不好，可以自己再調一調
beta = 3;              % weight of sparsity penalty term 稀疏自編碼稀疏項的權重

%%======================================================================
%% Load data from the MNIST database
%  This loads our training data from the MNIST database files.
% Load MNIST database files
mnistData = loadMNISTImages('train-images.idx3-ubyte');
mnistLabels = loadMNISTLabels('train-labels.idx1-ubyte');
mnistLabels(mnistLabels == 0) = 10; % Remap 0 to 10 since our labels need to start from 1
trainData   = mnistData;%訓練樣本
trainLabel= mnistLabels;

接下來我們先訓練第一個稀疏自編碼（這塊跟以前代碼一樣，就不注釋了，可以看以前博客的代碼，有注釋）：

%% Train the first sparse autoencoder
%  Randomly initialize the parameters
sae1Theta = initializeParameters(hiddenSizeL1, inputSize);
addpath minFunc/
options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost
                          % function. Generally, for minFunc to work, you
                          % need a function pointer with two outputs: the
                          % function value and the gradient. In our problem,
                          % sparseAutoencoderCost.m satisfies this.
options.maxIter = 400;    % Maximum number of iterations of L-BFGS to run
options.display = 'on';
[sae1OptTheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ...
                                   inputSize, hiddenSizeL1, ...
                                   lambda, sparsityParam, ...
                                   beta, trainData), ...
                              sae1Theta, options);
%現在訓練好了第一個的稀疏編碼神經網絡

我們從下圖可以看到最后cost function 收斂到 12.22附近

 

訓練好以后我們想直觀看看我們訓練出來的那200個特征（代碼為什么這么寫原因在以前博客中已經說明了）：

W1 = reshape(sae1OptTheta(1:hiddenSizeL1 * inputSize), hiddenSizeL1, inputSize);
display_network(W1');

稀疏自編碼學習的200個特征如下：

 

訓練第二個的稀疏編碼神經網絡：

%% Train the second sparse autoencoder
%開始訓練第二個的稀疏編碼神經網絡，注意輸入到第二個稀疏編碼神經網絡的訓練樣本應該是第一個稀疏編碼的激活值。
%所以要把第一個的訓練樣本unlabeledData輸入到feedForwardAutoencoder（）得到激活值，做為第二個稀疏編碼神經網絡的訓練樣本。
[sae1Features] = feedForwardAutoencoder(sae1OptTheta, hiddenSizeL1, ...
                                        inputSize, trainData);%sae1Features就是第二個稀疏編碼神經網絡的訓練樣本

%進行第二個的稀疏編碼神經網絡訓練，和第一個的方法一樣
sae2Theta = initializeParameters(hiddenSizeL2, hiddenSizeL1);%這時候訓練樣本是第一層隱藏層的激活值，所以
%這時候的'inputSize'是hiddenSizeL1。
addpath minFunc/
options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost
                          % function. Generally, for minFunc to work, you
                          % need a function pointer with two outputs: the
                          % function value and the gradient. In our problem,
                          % sparseAutoencoderCost.m satisfies this.
options.maxIter = 400;    % Maximum number of iterations of L-BFGS to run
options.display = 'on';
[sae2OptTheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ...
                                    hiddenSizeL1, hiddenSizeL2, ...
                                   lambda, sparsityParam, ...
                                   beta, sae1Features), ...
                              sae2Theta, options);
%現在訓練好了第二個的稀疏編碼神經網絡

我們從下圖可以看到最后cost function 收斂到 3.589附近

訓練好第二個的稀疏編碼神經網絡以后，我們想看下trainData的二階特征是什么樣的？為了說了更詳細一點，這個放到下一個博客當中。

接下來我們就要訓練整個模型的最后一層softmax 層了（這里再啰嗦幾句，雖然softmax和前面稀疏自編碼都需要+1的截距神經元，但是因為以前寫的函數里面后者包括了+1的稀疏自編碼，所以樣本有多少維、有多少神經元有多少參數，你就不用管+1截距神經元了。而softmax你需要管這些）：

%% Train the softmax classifier
%  This trains the sparse autoencoder on the second autoencoder features.
[sae2Features] = feedForwardAutoencoder(sae2OptTheta, hiddenSizeL2, ...
                                        hiddenSizeL1, sae1Features);%把sae1Features輸入到訓練好的第2個稀疏編碼
                                    %神經網絡中即使用feedForwardAutoencoder，得到sae2Features，也就是用來訓練
                                    %softmax模型的訓練樣本,sae2Features大小和以前的data數據格式一樣，每一列是一                  %個樣本。每一個樣本有hiddenSizeL2維。
sae2Features=[ones(1,size(sae2Features,2));sae2Features];%增加一維截距+1
%在這里你不需要初始化softmax模型參數，因為下面softmaxTrain函數里面會自動初始化，還有輸入softmaxTrain的‘putinsize’和樣本應該包括
%+1的截距神經元，所以這里是‘putinsize’是hiddenSizeL2+1
options.maxIter = 100;
softmaxModel = softmaxTrain(hiddenSizeL2+1, numClasses, lambda2, ...
                            sae2Features, trainLabel, options);
%注意softmaxTrain里面有softmax_regression_vec函數，softmax_regression_vec函數輸入的y需要是一個列向量（因為有sub2ind函數，具體看函  %數內部注釋），所以這里的trainlabel
%需要是列向量。
saeSoftmaxOptTheta = softmaxModel.optTheta(:);
%softmax模型的參數已經訓練好了。saeSoftmaxOptTheta是一個列向量。

我們從下圖可以看到最后cost function 收斂到 0.368附近

接下來就是重頭戲了，我們要微調整個網絡了：

我們先初始化微調這一步整個網絡的參數，注意這里不是隨機的給出參數了，因為是微調，所以初始化的值是之前訓練好的2個稀疏自編碼模型的參數和

softmax層的參數。

%% Finetune softmax model
% Initialize the stackedAETheta using the parameters learned
%我們從前面的博客中已經說了,我們訓練好每一個稀疏自編碼，我們需要的是每個神經網絡的第一層的參數，即w（1）和b（1），
%stack{1}是第一個稀疏自編碼里面的w（1）和b（1）參數,stack{2}標記的是第二個稀疏自編碼w（1）和b（1）參數
stack = cell(2,1);
stack{1}.w = reshape(sae1OptTheta(1:hiddenSizeL1*inputSize), ...
                     hiddenSizeL1, inputSize);
stack{1}.b = sae1OptTheta(2*hiddenSizeL1*inputSize+1:2*hiddenSizeL1*inputSize+hiddenSizeL1);
stack{2}.w = reshape(sae2OptTheta(1:hiddenSizeL2*hiddenSizeL1), ...
                     hiddenSizeL2, hiddenSizeL1);
stack{2}.b = sae2OptTheta(2*hiddenSizeL2*hiddenSizeL1+1:2*hiddenSizeL2*hiddenSizeL1+hiddenSizeL2);

% Initialize the parameters for the deep model
[stackparams, netconfig] = stack2params(stack);%把原胞stack變成一個向量，注意這個向量里面元素代表的是什么，順序不要混。
stackedAETheta = [ saeSoftmaxOptTheta ; stackparams ];%這個就是整個模型進行微調時參數的初始化
%stackedAETheta向量的頭(hiddenSizeL2 +1)* numClasses個參數 是softmax 模型的參數。
save('debug2') %把前面運行完的得到的所有變量save 一下，因為接下來我們要看看寫得wholeNetCost_and_grad函數 是否正確，
%又因為原模型太復雜，運行起來浪費時間，所以我們簡化了一下模型，因此模型的結構和參數值都變了。我們再檢測完正確以后，我們想訓練這個微調這一步。
%我們還是需要原來整個模型的值，所以我們在檢測這個函數正確以后，再load 一下debug2.

上面程序有一個stack2params函數：

function [params, netconfig] = stack2params(stack)

% Converts a "stack" structure into a flattened parameter vector and also
% stores the network configuration. This is useful when working with
% optimization toolboxes such as minFunc.
%
% [params, netconfig] = stack2params(stack)
%
% stack - the stack structure, where stack{1}.w = weights of first layer
%                                    stack{1}.b = weights of first layer
%                                    stack{2}.w = weights of second layer
%                                    stack{2}.b = weights of second layer
%                                    ... etc.


% Setup the compressed param vector
params = [];
for d = 1:numel(stack)

    % This can be optimized. But since our stacks are relatively short, it
    % is okay
    params = [params ; stack{d}.w(:) ; stack{d}.b(:) ];

    % Check that stack is of the correct form
    assert(size(stack{d}.w, 1) == size(stack{d}.b, 1), ... %assert（expression,'string'）若expression是false的，那么就
        ['The bias should be a *column* vector of ' ...                              報錯string那段話
         int2str(size(stack{d}.w, 1)) 'x1']);
    if d < numel(stack)
        assert(size(stack{d}.w, 1) == size(stack{d+1}.w, 2), ...
            ['The adjacent layers L' int2str(d) ' and L' int2str(d+1) ...
             ' should have matching sizes.']);
    end

end
%params最后得到的是一個列向量，列向量的前hiddenSizeL1*inputSize個是以前訓練的第一個稀疏自編碼的w（1）（這里不包括b（1）），
%也是新的'大神經網絡'(上面的圖2)的第一層的參數（不包括+1的截距參數），列向量的第hiddenSizeL1*inputSize+1到
%hiddenSizeL1*inputSize+hiddenSizeL1是以前訓練的第一個稀疏自編碼的b（1）。列向量的第
%hiddenSizeL1*inputSize+hiddenSizeL1+1到hiddenSizeL1*inputSize+hiddenSizeL1+hiddenSizeL2*hiddenSizeL1
%以前訓練的第二個稀疏自編碼的w（1）（這里不包括b（1）），也是新的'大神經網絡'(上面的圖2)的第二層的參數（不包括+1的截距參數）
%第hiddenSizeL1*inputSize+hiddenSizeL1+hiddenSizeL2*hiddenSizeL1+1
%到hiddenSizeL1*inputSize+hiddenSizeL1+hiddenSizeL2*hiddenSizeL1+hiddenSizeL2是以前訓練的第二個稀疏自編碼的b（1）。
%但是我們發現這個‘大神經網絡’第一、二層參數w和截距參數b在params的排列跟以前的神經網絡的參數向量排列不怎么一樣。(以前的把所有
%層的+1的截止參數都放最后)這個注意一下。
if nargout > 1
    % Setup netconfig
    if numel(stack) == 0
        netconfig.inputsize = 0;
        netconfig.layersizes = {};
    else
        netconfig.inputsize = size(stack{1}.w, 2);
        netconfig.layersizes = {};
        for d = 1:numel(stack)
            netconfig.layersizes = [netconfig.layersizes ; size(stack{d}.w,1)];
        end
    end
end

end

下面看看寫得wholeNetCost_and_grad函數，這個函數的作用就是給定模型的參數θ，把樣本都通過這個總的模型得到相應的模型輸出，根據樣本的label這個函數就給出了當前模型參數下的 cost function 和各模型參數的導數grad。

function [cost,grad] = wholeNetCost_and_grad(theta, visibleSize,hiddenSizeL1,hiddenSizeL2,numClasses ,...
                                             lambda,lambda2, sparsityParam, beta, data,label)
%這個wholeNetCost_and_grad函數的作用是：你給定整個神經網絡的結構，網絡參數值，L2正則項、稀疏約束項的權重lambda，
%softmax的正則項lambda2，稀疏性參數。（在結構中注意一點，因為hiddenSizeL2即作為前面神經網絡的輸出，又作為最后
%一層softmax的輸入，而做為前面神經網絡的輸出時，用的是以前寫得稀疏自編碼的程序（都是三層，只是第三層神經元個數不再是稀疏自編碼
%輸入神經元的個數，而是hiddenSizeL2），所以考慮每一層神經元個數或參數的時候，是不用考慮+1的截距神經元的，不是說稀疏自編碼不需要，
%而是以前的函數內部已經有了，你只要輸入除+1截距神經元的個數就可以了，而對於softmax層的輸入，我們考慮ipuSixe的時候或有多少參數的
%時候，你要加上+1的截距神經元。所以下面對於softmax 的輸入我們要hiddenSizeL2+1）
%那么這個函數就會給出這個神經網絡的cost function的值（整個網絡的cost function），和cost function對各個參數的偏導。
%注意label輸入需要是一個列向量
%下面的W1，W2，b1，b2是總的模型第一層第二層參數，注意他在theta中的位置是什么，要細心。
W1 = reshape(theta((hiddenSizeL2 +1)* numClasses+1:(hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize),...
                                                                            hiddenSizeL1, visibleSize);
W2 = reshape(theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+1:...
(hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+ hiddenSizeL2*hiddenSizeL1),...
                                                                         hiddenSizeL2, hiddenSizeL1);
b1 = theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+1:...
                         (hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1);
b2 = theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+ hiddenSizeL2*hiddenSizeL1+1:end);
theta_softmax=reshape(theta(1:(hiddenSizeL2 +1)* numClasses),(hiddenSizeL2 +1),numClasses);
theta_softmax=theta_softmax(2:end,:);%求倒數第二層神經元殘差時要用到最后一層的參數θ矩陣，但是這里面的θ矩陣要去掉截距參數
%應該倒數第二層+1的神經元不用求殘差。
%% 初始化要求得的cost function 對各參數的偏導數
% Cost and gradient variables (your code needs to compute these values).
% Here, we initialize them to zeros.
cost = 0;
W1grad = zeros(size(W1)); %W1grad 也是一個矩陣，跟W1的行列數是一樣的。里面的元素與神經網絡里面的參數對應的情況和W1
%里元素與神經網絡參數對應的情況是一樣的。第一列各元素的含義是cost function 對輸入層第一個神經元分別對應第二層隱藏層的
%各單元（不包括第二層偏置節點）的參數的偏導。
W2grad = zeros(size(W2));
b1grad = zeros(size(b1));
b2grad = zeros(size(b2));

%% 求cost function 對各參數的偏導數
%  Instructions: Compute the cost/optimization objective J_sparse(W,b) for the Sparse Autoencoder,
%                and the corresponding gradients W1grad, W2grad, b1grad, b2grad.
%
% W1grad, W2grad, b1grad and b2grad should be computed using backpropagation.
% Note that W1grad has the same dimensions as W1, b1grad has the same dimensions
% as b1, etc.  Your code should set W1grad to be the partial derivative of J_sparse(W,b) with
% respect to W1.  I.e., W1grad(i,j) should be the partial derivative of J_sparse(W,b)
% with respect to the input parameter W1(i,j).  Thus, W1grad should be equal to the term
% [(1/m) \Delta W^{(1)} + \lambda W^{(1)}] in the last block of pseudo-code in Section 2.2
% of the lecture notes (and similarly for W2grad, b1grad, b2grad).
%
% Stated differently, if we were using batch gradient descent to optimize the parameters,
% the gradient descent update to W1 would be W1 := W1 - alpha * W1grad, and similarly for W2, b1, b2.
%
%1.forward propagation
data_size=size(data);
biasPara_1=repmat(b1,1,data_size(2));
biasPara_2=repmat(b2,1,data_size(2));
active_value2=sigmoid(W1*data+biasPara_1);
active_value3=sigmoid(W2*active_value2+biasPara_2);
%active_value2這個矩陣每一列是每一個樣本前向算法中第二層各單元的激活值（不包括偏置節點，因為偏置節點
%相當於激活值始終為1）。
%active_value3是第三層各單元的激活值。

%最后一層的'傳播'：把active_value3增廣數據輸入softmax中得到這個模型最后的輸出，使用的是softmax_regression_vec_addP函數
%softmax_regression_vec_addP是softmax_regression_vec在輸出的時候增加了一個p，p矩陣每一列是每個樣本對應的各個輸出神經元的概率
softmax_input=[ones(1,size(active_value3,2));active_value3];%每個樣本都加一個+1維
[p,cost,softmax_grad] = softmax_regression_vec_addP(theta(1:(hiddenSizeL2 +1)* numClasses),softmax_input ,...
                                                                                           label,lambda2 );
%這里的softmax_grad是一個向量，是softmax這一層的參數的導數。
weight_decay=lambda/2*(sum(sum(W1.^2))+sum(sum(W2.^2)));
%這是const function 中前面幾層網絡的L2正則項，softmax
%的L2正則項在softmax_regression_vec_addP已經有了，所以上面的那個cost 包括了。
active_average=[sum(active_value2,2)./data_size(2);sum(active_value3,2)./data_size(2)];
%因為這里有兩層隱藏層，並且active_average是一個列向量的形式，里面每一個元素是每一個隱藏層神經元在所有樣本下的平均激活值
p_para=repmat(sparsityParam,hiddenSizeL1+hiddenSizeL2,1);
%因為下面那行代碼中有p_para./active_average所以要弄成和active_average一樣size的列向量，才能與active_average里的元素點除。
sparsity=beta.*sum(p_para.*log(p_para./active_average)+(1-p_para).*log((1-p_para)./(1-active_average)));
%這是稀疏項
cost=cost+weight_decay+sparsity;%cost function就完成了。
%現在我們有了整個模型在特定θ下的cost function，也得到了最后一層softmax參數的導數，只剩下前面幾層神經網絡的導數了，
%現在我們用后向傳播算法和一些其他公式進行求出，這里跟以前稍有不同的是整個模型倒數第二層神經元的殘差，不同的地方已經在上面博客中
%說明了。

%求殘差
%根據博客中公式（稍有不同的是，這里是所有樣本一塊操作，每一列是一個樣本的信息）
I=full(sparse(label,1:data_size(2),1)); %I為一個稀疏矩陣，這個矩陣的（label（1）,1）、（label（2）,2）
  %....(label(m),m)位置都是1。
  active_average_repmat1=repmat(sum(active_value3,2)./data_size(2),1,data_size(2));
  default_sparsity1=repmat(sparsityParam,hiddenSizeL2,data_size(2));
  sparsity_penalty1=beta.*(-(default_sparsity1./active_average_repmat1)+((1-default_sparsity1)./...
                                                             (1-active_average_repmat1)));
delta3=(theta_softmax*(p-I)+sparsity_penalty1).*(active_value3).*(1-active_value3);
%因為這里delta3、theta_softmax*(p-I)、active_value3都是一個矩陣，每一列是一個樣本輸入神經網絡后第三層各個神經元的殘差、激活值、小cost
%function（單樣本）對第三層各個神經元的導數，列數代表樣本數。所以sparsity_penalty1、active_average_repmat1、sparsity_penalty1
%要和前面提的矩陣大小相等。所以default_sparsity1是平鋪hiddenSizeL2行。下面的default_sparsity也是同樣道理平鋪hiddenSizeL1行
active_average_repmat=repmat(sum(active_value2,2)./data_size(2),1,data_size(2));
%sum(active_value2,2)./data_size(2)就是active_average，然后repmat，因為active_average是一個列向量，再平鋪
%data_size(2)列，那么active_average_repmat就是一個矩陣，弄成矩陣的的形式是為了后面所有樣本第二層所有神經元
%一起運算。
default_sparsity=repmat(sparsityParam,hiddenSizeL1,data_size(2));
%default_sparsity也是一個hiddenSize*data_size(2)矩陣，弄成這個形式也是為了下面所有樣本第二層所有神經元一起的矩陣
%運算。
sparsity_penalty=beta.*(-(default_sparsity./active_average_repmat)+((1-default_sparsity)./(1-active_average_repmat)));
delta2=(W2'*delta3+sparsity_penalty).*((active_value2).*(1-active_value2));

%求各參數偏導
W2grad=delta3*active_value2'./data_size(2)+lambda.*W2;
%delta3*active_value2'得到的是一個矩陣，矩陣的第一個元素是所有樣本對應的 小cost function
%對第二層的W11的偏導 的累加形成的總cost function對第二層的W11的偏導。矩陣第一行第二個元素是
%總cost function對第二層的W21的偏導。所以你看這里的W下標的順序和前面的W矩陣是一樣的。然后這個矩陣的元素都除以樣本
%個數，得到的就是cost function中只有均方誤差項時，對第二層各個參數的偏導。然后這個矩陣每個元素加上lambda倍的W2這個
%矩陣對應的元素，得到的W2grad矩陣的元素是總cost function對第二層各參數的偏導。
W1grad=delta2*data'./data_size(2)+lambda.*W1;
b2grad=sum(delta3,2)./data_size(2);%b2grad向量是總cost function對第二層偏置節點各參數的偏導
b1grad=sum(delta2,2)./data_size(2);

%-------------------------------------------------------------------
% After computing the cost and gradient, we will convert the gradients back
% to a vector format (suitable for minFunc).  Specifically, we will unroll
% your gradient matrices into a vector.

grad = [softmax_grad;W1grad(:) ;b1grad(:) ; W2grad(:) ;  b2grad(:)];%注意里面元素的順序和theta是一樣的。

end

%-------------------------------------------------------------------
% Here's an implementation of the sigmoid function, which you may find useful
% in your computation of the costs and the gradients.  This inputs a (row or
% column) vector (say (z1, z2, z3)) and returns (f(z1), f(z2), f(z3)).
%若果x是個矩陣，也是可以的，igmoid(x)也是一個矩陣，里面的每個元素就是對應x矩陣對應的元素進行sigmoid計算。

function sigm = sigmoid(x)

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

接下面我們要對寫的這個函數檢驗一下，在給定模型參數下，模型的cost function 對各參數求導是不是正確，檢驗方法是利用求導的定義，在前面已經講了，因為這個模型和參數比較多，我們簡化一下模型，對檢驗這個函數正確性上沒有影響。

%% 看看寫得wholeNetCost_and_grad 是否正確
DEBUG = true;
if DEBUG
    stackedAETheta=0.005 * randn(430, 1);
    inputSize = 20;
    trainData = randn(20, 100);
    hiddenSizeL1=10;
    hiddenSizeL2=10;
    numClasses=10;
    trainLabel= randi(10, 100, 1);%從[1,10]中隨機生成一個100*1的列向量
[cost,grad] = wholeNetCost_and_grad(stackedAETheta, inputSize, hiddenSizeL1,hiddenSizeL2,numClasses, ...
                                             lambda,lambda2, sparsityParam, beta,trainData,trainLabel);
     %注意trainlabel輸入需要是一個列向量


numGrad = computeNumericalGradient( @(x) wholeNetCost_and_grad(x, inputSize, hiddenSizeL1,...
                                                            hiddenSizeL2,numClasses, ...
                                             lambda,lambda2, sparsityParam, beta,trainData,trainLabel),...
                                       stackedAETheta  );

disp([numGrad grad]);
diff = norm(numGrad-grad)/norm(numGrad+grad);
disp(diff);
end

結果為：

所以我們知道我們寫的函數是正確的。

那么接下來就是微調了：

%% 微調 訓練 前面的stackedAETheta就是這里最小化過程中的初始化參數向量
load debug2  %把前面的變量值再加載一遍
addpath minFunc/
options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost
                          % function. Generally, for minFunc to work, you
                          % need a function pointer with two outputs: the
                          % function value and the gradient. In our problem,
                          % sparseAutoencoderCost.m satisfies this.
options.maxIter = 400;    % Maximum number of iterations of L-BFGS to run
options.display = 'on';
[stackedAEOptTheta, cost] = minFunc( @(p) wholeNetCost_and_grad(p, ...
                                     inputSize, hiddenSizeL1,hiddenSizeL2,numClasses, ...
                                             lambda,lambda2, sparsityParam, beta,trainData,trainLabel), ...
                              stackedAETheta, options);
%stackedAEOptTheta就是最終整個模型的參數。

現在模型的參數都已經訓練好了，接下來就是檢驗了（檢驗用的數據庫前面博客中已經給出了，可以去那里下載）：

%% 測試
testData = loadMNISTImages('t10k-images.idx3-ubyte');
testLabels = loadMNISTLabels('t10k-labels.idx1-ubyte');

testLabels(testLabels == 0) = 10; % Remap 0 to 10

[pred] = stackedAEPredict(stackedAETheta, inputSize, hiddenSizeL1,hiddenSizeL2, numClasses,  testData);

acc = mean(testLabels(:) == pred(:));
fprintf('Before Finetuning Test Accuracy: %0.3f%%\n', acc * 100);

[pred] = stackedAEPredict(stackedAEOptTheta, inputSize,  hiddenSizeL1,hiddenSizeL2, numClasses,  testData);

acc = mean(testLabels(:) == pred(:));
fprintf('After Finetuning Test Accuracy: %0.3f%%\n', acc * 100);

代碼中的stackedAEPredict函數：

function [pred] = stackedAEPredict(theta, visibleSize, hiddenSizeL1,hiddenSizeL2, numClasses, data)

% stackedAEPredict: Takes a trained theta and a test data set,
% and returns the predicted labels for each example.

% theta: trained weights from the autoencoder
% visibleSize: the number of input units
% hiddenSize:  the number of hidden units *at the 2nd layer*
% numClasses:  the number of categories
% data: Our matrix containing the training data as columns.  So, data(:,i) is the i-th training example.

% Your code should produce the prediction matrix
% pred, where pred(i) is argmax_c P(y(c) | x(i)).

%% 下面的代碼都是利用上面函數的部分代碼，道理都是一樣的，不懂的可以往前看看
W1 = reshape(theta((hiddenSizeL2 +1)* numClasses+1:(hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize),...
                                                                            hiddenSizeL1, visibleSize);
W2 = reshape(theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+1:...
(hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+ hiddenSizeL2*hiddenSizeL1),...
                                                                         hiddenSizeL2, hiddenSizeL1);
b1 = theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+1:...
                         (hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1);
b2 = theta((hiddenSizeL2 +1)* numClasses+hiddenSizeL1*visibleSize+hiddenSizeL1+ hiddenSizeL2*hiddenSizeL1+1:end);
data_size=size(data);
biasPara_1=repmat(b1,1,data_size(2));
biasPara_2=repmat(b2,1,data_size(2));
active_value2=sigmoid(W1*data+biasPara_1);
active_value3=sigmoid(W2*active_value2+biasPara_2);
softmax_input=[ones(1,size(active_value3,2));active_value3];%每個樣本都加一個+1維
p= softmax_regression_vec_justP(theta(1:(hiddenSizeL2 +1)* numClasses),softmax_input );%其實用的還是softmax_regression_vec
%函數，只是這里讓它僅輸出p
[q,pred]=max(p);%這里認為最大概率的類別就是樣本類別
end

% You might find this useful
function sigm = sigmoid(x)
    sigm = 1 ./ (1 + exp(-x));
end

我們從下圖可以看到最后cost function 收斂到 0.655附近

全部就完成了，看下我們最終的結果吧：

我們發現微調以后我們的正確率提高了3%多一點，我們可以通過調參數進一步的提高正確率。