本文分為三部分，前言和實驗基礎copy自tornadomeet部落格；後面代碼詳解為自己原創。
前言：
本次是練習2個隐含層的網絡的訓練方法，每個網絡層都是用的sparse autoencoder思想，利用兩個隐含層的網絡來提取出輸入資料的特征。本次實驗驗要完成的任務是對MINST進行手寫數字識别，實驗内容及步驟參考網頁教程Exercise: Implement deep networks for digit classification。當提取出手寫數字圖檔的特征後，就用softmax進行對其進行分類。關于MINST的介紹可以參考網頁：MNIST Dataset。本文的理論介紹也可以參考前面的博文：Deep learning：十六(deep networks)。
實驗基礎：
進行deep network的訓練方法大緻如下：
1. 用原始輸入資料作為輸入，訓練出（利用sparse autoencoder方法）第一個隐含層結構的網絡參數，并将用訓練好的參數算出第1個隐含層的輸出。
2. 把步驟1的輸出作為第2個網絡的輸入，用同樣的方法訓練第2個隐含層網絡的參數。
3. 用步驟2 的輸出作為多分類器softmax的輸入，然後利用原始資料的标簽來訓練出softmax分類器的網絡參數。
4. 計算2個隐含層加softmax分類器整個網絡一起的損失函數，以及整個網絡對每個參數的偏導函數值。
5. 用步驟1，2和3的網絡參數作為整個深度網絡（2個隐含層,1個softmax輸出層）參數初始化的值，然後用lbfs算法疊代求出上面損失函數最小值附近處的參數值，并作為整個網絡最後的最優參數值。
上面的訓練過程是針對使用softmax分類器進行的，而softmax分類器的損失函數等是有公式進行計算的。是以在進行參數校正時，可以對把所有網絡看做是一個整體，然後計算整個網絡的損失函數和其偏導，這樣的話當我們有了标注好了的資料後，就可以用前面訓練好了的參數作為初始參數，然後用優化算法求得整個網絡的參數了。但如果我們後面的分類器不是用的softmax分類器，而是用的其它的，比如svm，随機森林等，這個時候前面特征提取的網絡參數已經預訓練好了，用該參數是可以初始化前面的網絡，但是此時該怎麼微調呢？因為此時标注的數值隻能在後面的分類器中才用得到，是以沒法計算系統的損失函數等。難道又要将前面n層網絡的最終輸出等價于第一層網絡的輸入（也就是多網絡的sparse autoencoder）?本人暫時還沒弄清楚，日後應該會想明白的。
關于深度網絡的學習幾個需要注意的小點（假設隐含層為2層）：
- 利用sparse autoencoder進行預訓練時，需要依次計算出每個隐含層的輸出，如果後面是采用softmax分類器的話，則同樣也需要用最後一個隐含層的輸出作為softmax的輸入來訓練softmax的網絡參數。
- 由步驟1可知，在進行參數校正之前是需要對分類器的參數進行預訓練的。且在進行參數校正(Finetuning )時是将所有的隐含層看做是一個單一的網絡層，是以每一次疊代就可以更新所有網絡層的參數。
另外在實際的訓練過程中可以看到，訓練第一個隐含層所用的時間較長，應該需要訓練的參數矩陣為200*784(沒包括b參數),訓練第二個隐含層的時間較第一個隐含層要短些，主要原因是此時隻需學習到200*200的參數矩陣，其參數個數大大減小。而訓練softmax的時間更短，那是因為它的參數個數更少，且損失函數和偏導的計算公式也沒有前面兩層的複雜。最後對整個網絡的微調所用的時間和第二個隐含層的訓練時間長短差不多。
以上複制來自：http://www.cnblogs.com/tornadomeet/archive/2013/04/09/3011209.html
代碼流程詳解：
Step0:
初始化網絡參數，确定網絡結構
28*28 100 100 10
Inputsize hiddensizeL1 hiddensizeL1 numclass
W1=matrix（hiddensizeL1，inputsize）= matrix（100,764）
b1=vector（hiddensizeL1，1）=vector（100,1）
sae1Theta=[W1,b1];
W2=matrix（hiddensizeL2，hiddensizeL1）= matrix（100,100）
b2=vector（hiddensizeL2，1）=vector（100,1）
SoftmaxTheta=matrix（numclass，hiddensizeL2）=matrix（10,100）
Step1:加載資料
trainData = loadMNISTImages('../data/train-images-idx3-ubyte');
Step2：訓練第一個稀疏自編碼，得到W1和b1的優化權值
[sae1OptTheta,loss] = minFunc( @(p) sparseAutoencoderCost(p, ...
利用函數sparseAutoencoderCost（）計算損失值和梯度；并利用minFunc的l-bfgs
算法進行優化，得到優化後的W1和b1的權值，sae1OptTheta=[W1,b1]。
Step3：訓練第二個稀疏自編碼，得到W2和b2的優化權值
1.
[sae1Features] = feedForwardAutoencoder(sae1OptTheta, ..
利用前饋函數feedForwardAutoencoder計算原始資料data的第一個特征表達sae1Features
2.
初始化W2和b2
sae2Theta = initializeParameters(hiddenSizeL2,hiddenSizeL1);
3.
利用函數sparseAutoencoderCost（）計算損失值和梯度；并利用minFunc的算法進行優化，得到優化後的W2和b2的權值，sae2OptTheta=[W2,b2]。
Step4:訓練softmax分類器
1.
[sae2Features] = feedForwardAutoencoder(sae2OptTheta,hiddenSizeL2, ...hiddenSizeL1, sae1Features);
利用前饋函數和第一次特征表達sae1Features計算第二層特征表達sae2Features 作為softmax分類器的輸入資料
2.
初始化softmax分類器的參數theta，用小數值随機初始化
saeSoftmaxTheta = 0.005 * randn(hiddenSizeL2 *numClasses, 1);
3.
訓練softmax分類器
softmaxModel = softmaxTrain(hiddenSizeL2, numClasses,1e-4, ... sae2Features,trainLabels, options);
softmaxTrain函數的輸出參數softmaxModel為一個結構體，即struct；其中包含三個子屬性值
softmaxModel.optTheta =reshape(softmaxOptTheta, numClasses, inputSize); %分類器的參數theta
softmaxModel.inputSize = inputSize; %分類器的輸入參數個數
softmaxModel.numClasses = numClasses; %分類器的種類
3.1
優化參數theta
在softtrain中調用的函數softmaxcost，類似sparseAutoencoderCost函數，通過softmax分類器模型來計算函數損失值和梯度值，之後通過l-bfgs算法來優化。
Step5：微調fine-tuning softmax模型
Setp5就是這個算法的關鍵部分了，首先這個算法叫做stackedautoencoder，棧式自編碼，其中這個“棧”，也就是stack，就展現在這裡。這裡的stack就是一個元包（cell），cell的元素是各個參數矩陣，即W1，b1等。通過stack的初始化，把這個深度網絡結構的所有參數W1，b1，W2，b2合并（打包）到stack中；後續計算中通過打包（params2stack函數完成）與拆包（stack2params）函數完成)來計算更新各個權值。
5.1 參數打包到stack中，進行棧化
stack = cell(2,1);%定義一個包含2個元包的元包矩陣（額，稱呼有些不恰當）
stack{1}.w =reshape(sae1OptTheta(1:hiddenSizeL1*inputSize), ...hiddenSizeL1, inputSize);
stack{1}，是調用stack中的第一個元包，第一個元包的資料類型是結構體struct；裡面包含W1和b1兩個屬性值。
stack{1}.w = weights of first layer=W1
stack{1}.b = weights of first layer=b1
stack{2}.w = weights of second layer=W2
stack{2}.b = weights of second layer=b2
[stackparams, netconfig] = stack2params(stack);
stackedAETheta = [ saeSoftmaxOptTheta ; stackparams ];%合并網絡參數和分類器參數。
5.1.1 stack2params函數說明
[stackparams, netconfig] =stack2params(stack);
把棧stack轉換為params的過程，其輸出參數為兩個，一個是向量化後的網絡參數stackparams，一個是網絡配置結構參數netconfig；這種方式感覺很簡潔，應為我們訓練神經網絡最後要得到的就是網絡結構和結構件的連接配接權重，stackparams來記錄網絡的連接配接權重，netconfig來記錄網絡的配置結構；這樣就把一個複雜的網絡，通過兩個略微複雜的變量記錄下來，友善後續的應用和計算。
函數部分一：
打包過程，stack向量化
for d = 1:numel(stack)
params = [params ; stack{d}.w(:) ; stack{d}.b(:) ];
end
assert函數說明：
MATLAB語言沒有系統的斷言函數，但有錯誤報告函數error 和 warning。由于要求對參數的保護，需要對輸入參數或處理過程中的一些狀态進行判斷，判斷程式能否/是否需要繼續執行。
其中的assert函數，類似if函數，判斷程式是否進行
assert(size(stack{d}.w, 1) == size(stack{d}.b, 1), ['The biasshould be a *column* vector of ' int2str(size(stack{d}.w,1)) 'x1']);
等價于：
If size(stack{d}.w,1) ~= size(stack{d}.b, 1)
temp= int2str(size(stack{d}.w, 1))；
error（‘The biasshould be a *column* vector of ‘temp’x1’）
end
用來判斷stack{d}.w的行數和stack{d}.b的行數是否相同，都等于隐層節點數目。
函數部分二：
netconfig（netconfiguration，網絡配置）部分，記錄網絡各層的配置參數（節點數）。
變量nargout，是matlab函數自帶的“隐藏”變量，判斷函數中輸出參數個數的變量，n arg out
if nargout> 1 …… end %如果本函數的輸出參數大于一個，則執行下面的代碼
netconfig變量是一個struct變量，有inputsize（一個數值）和layersizes（一個量）向兩個子屬性。
netconfig.inputsize = size(stack{1}.w, 2);netconfig.layersizes= {};
for d = 1:numel(stack)
netconfig.layersizes = [netconfig.layersizes ;size(stack{d}.w,1)];
end
周遊stack的所有元包，通過size(stack{d}.w,1)，來記錄網絡每層的節點數。
5.2
利用stackedAECost函數計算梯度和函數損失值，并用l-bfgs算法優化
stackedAECost函數說明
1.把向量化參數theta還原為矩陣形式
softmaxTheta =reshape(theta(1:hiddenSize*numClasses), numClasses, hiddenSize);
stack =params2stack(theta(hiddenSize*numClasses+1:end), netconfig);
2.初始化梯度矩陣softmaxThetaGrad和梯度棧stackgrad
softmaxThetaGrad = zeros(size(softmaxTheta));
stackgrad = cell(size(stack));
for d = 1:numel(stack)
stackgrad{d}.w = zeros(size(stack{d}.w));
stackgrad{d}.b = zeros(size(stack{d}.b));
end
3.前饋計算神經網絡
depth = numel(stack);%擷取棧中元包個數，即網絡的層數
z = cell(depth+1,1);
a = cell(depth+1, 1);%記錄原始輸入資料和各層的輸出激活值a
a{1} = data;
for layer = (1:depth)
z{layer+1} = stack{layer}.w * a{layer} + repmat(stack{layer}.b, [1,size(a{layer},2)]);%前饋計算w*x+b
a{layer+1} = sigmoid(z{layer+1});%利用sigmoid函數進行非線性變換
end
4.計算softmax分類器
M = softmaxTheta * a{depth+1}; %計算theta*x
M = bsxfun(@minus, M, max(M));%參數theta備援消除，就是把參數變小
p = bsxfun(@rdivide, exp(M), sum(exp(M)));%機率歸一化操作
cost = -1/numCases * groundTruth(:)' *log(p(:)) + lambda/2 * sum(softmaxTheta(:) .^ 2);%計算softmax函數值y
softmaxThetaGrad = -1/numCases * (groundTruth- p) * a{depth+1}' + lambda * softmaxTheta;%計算梯度值，并更新
5.反向傳播，計算各層梯度
d = cell(depth+1);%存放梯度的元包
d{depth+1} = -(softmaxTheta' * (groundTruth -p)) .* a{depth+1} .* (1-a{depth+1});%計算最後一層（softmax）的梯度值
for layer = (depth:-1:2)
d{layer} = (stack{layer}.w' * d{layer+1}) .* a{layer} .* (1-a{layer});%通過反向傳播機制計算各層的梯度
end
6.梯度跟新
for layer = (depth:-1:1)
stackgrad{layer}.w = (1/numCases) * d{layer+1} * a{layer}';%梯度更新
stackgrad{layer}.b = (1/numCases) * sum(d{layer+1}, 2);
end
Step6：預測精确度
1.加載測試資料
testData = loadMNISTImages('../data/t10k-images-idx3-ubyte');
2.用沒有經過微調的參數進行預測
[pred] = stackedAEPredict(stackedAETheta, inputSize, hiddenSizeL2, ...numClasses,netconfig, testData);
stackedAEPredict（）函數說明
部分一：拆包，把參數矩陣化
softmaxTheta = reshape(theta(1:hiddenSize*numClasses), numClasses,hiddenSize);
stack= params2stack(theta(hiddenSize*numClasses+1:end), netconfig);
部分二：前饋網絡計算
for layer = (1:depth)
z{layer+1} = stack{layer}.w * a{layer} +repmat(stack{layer}.b, [1, size(a{layer},2)]);
a{layer+1} = sigmoid(z{layer+1});
end
[~,pred] = max(softmaxTheta * a{depth+1});
%pred=max（a），結果pred是矩陣a每列的最大值；[~, pred]=max（a），則pred是矩陣a每列最大值的位置索引，此處為對應的類型
3.預測輸出
acc = mean(testLabels(:) ==pred(:));
fprintf('Before Finetuning Test Accuracy: %0.3f%%\n', acc * 100);