UFLDL深度學習筆記 （一）反向傳播與稀疏自編碼

UFLDL深度學習筆記 （一）基本知識與稀疏自編碼

前言

　　近來正在系統研究一下深度學習，作為新入門者，為了更好地理解、交流，准備把學習過程總結記錄下來。最開始的規划是先學習理論推導；然后學習一兩種開源框架；第三是進階調優、加速技巧。越往后越要帶着工作中的實際問題去做，而不能是空中樓閣式沉迷在理論資料的舊數據中。深度學習領域大牛吳恩達(Andrew Ng)老師的UFLDL教程 (Unsupervised Feature Learning and Deep Learning)提供了很好的基礎理論推導，文辭既系統完備，又明白曉暢，最可貴的是能把在斯坦福大學中的教學推廣給全世界，盡心力而廣育人，實教授之典范。
　　這里的學習筆記不是為了重復UFLDL中的推導過程，而是為了一方面補充我沒有很快看明白的理論要點，另一方面基於我對課后練習的matlab實現，記錄討論卡殼時間較長的地方。也方便其他學習者，學習過程中UFLDL所有留白的代碼模塊全是自己編寫的，沒有查網上的代碼，經過一番調試結果都達到了練習給出的參考標准，而且都按照矩陣化實現(UFLDL稱為矢量化 vectorization)，代碼鏈接見此https://github.com/codgeek/deeplearning，所以各種matlab實現細節、緣由也會比較清楚，此外從實現每次練習的代碼commit中可以清楚看到修改了那些代碼，歡迎大家一起討論共同進步。

理論推導中的要點

先上一個經典的自編碼神經網絡結構圖。推導中常用的符號表征見於此

上圖的神經網絡結構中，輸入單個向量$\ x\ \(,逐層根據網絡參數矩陣\)\ W, b\ \(計算[**前向傳播**](http://deeplearning.stanford.edu/wiki/index.php/%E7%A5%9E%E7%BB%8F%E7%BD%91%E7%BB%9C)，最后輸出向量\) h_{W,b}(x)$, 這樣單個向量的誤差代價函數是

\[J{(W,b,x,y)}=\frac{1}{2}{||(h_{W,b}(x)-y||}^2 \]

對所有數據代價之和即

\[\begin{align} J(W,b) &= \frac 1 m \left [\sum_{i=1}^m J{(W,b;x^{(i)},y^{(i)})} \right] +\frac \lambda 2 \sum_{l=1}^{n_l-1}\sum_{i=1}^{s_l}\sum_{j=1}^{s_{l+1}}{W_{ij}^{(l)}}^2 \\ &= \frac 1 m \left [\sum_{i=1}^m{\frac 1 2||(h_{W,b}(x^{(i)})-y^{(i)}||}^2 \right] +\frac \lambda 2 \sum_{l=1}^{n_l-1}\sum_{i=1}^{s_l}\sum_{j=1}^{s_{l+1}}{W_{ij}^{(l)}}^2 \end{align} \]

偏導與反向傳導的前因后果

神經網絡的最基礎的理論是反向傳導算法，反向傳導的主語是誤差代價函數$\ J_{W,b}\ $對網絡系數(W，b)的偏導數。那為什么要求偏導呢? 這里就引出了最優化理論，具體的問題一般被建模成許多自變量決定的代價函數，最優化的目標就是找到使代價函數最小化的自變量取值。數學方法告訴我們令偏導為0對應的自變量取值就是最小化代價的點，不幸的是多數實際問題對用偏導為0方程求不出閉合的公式解，神經網絡也不例外。為了解決這類問題，可以使用優化領域的梯度下降法，梯度是對每個自變量的偏導組成的向量或矩陣，梯度下降直白地說就是沿着對每個自變量偏導的負方向改變自變量，下降指的是負方向，直到梯度值接近0或代價函數變化很小，即停止查找認為找到了達到最優化的所有自變量數值。 所以為了最小化稀疏自編碼誤差代價函數，我們需要求出J(W,b)對神經網絡各層系數W、b的偏導數組成的梯度矩陣。

接下來在反向傳導算法文章中給出了詳細的推導過程，推導的結果便是殘差、梯度值從輸出層到隱藏層的反向傳導公式

輸出層的對輸出值$\ z\ \(第i個分量\)\ z_i^{n_l}\ $的偏導為

\[\delta_i^{(n_l)}= \frac {\partial} {\partial z_i^{(n_l)}} {\frac 1 2||(h_{W,b}(x^{(i)})-y^{(i)}||}^2 = -(y_i - a_i^{(n_l)}) \cdot f'( z_i^{(n_l)}) \]

后向傳導時，第l+1層傳導給第l層偏導值為 $$\delta_i^{(l)} = (\sum_{j=1}^{s_l+1}W_{ij}^{(l)}\delta_j^{(l+1)}) f'( z_i^{(l)}) $$

J對輸出值z的偏導怎么得到對網絡權值參數的偏導？

我們知道$\ z_i^{l+1}\ \(與\)\ W,b\ \(的關系\)\ z_i^{(l+1)}=\sum_{j=1}^{n_l} ({W_{i,j}^{(l)}\cdot{a_j^{(l)}+b}})\ $，運用求導鏈式法則可得：

\[\begin{align} \frac {\partial J_{W,b;x,y}}{\partial W_{i,j}^{(l)}} &= \frac {\partial J_{W,b;x,y}}{\partial z_i^{(l+1)}} \cdot \frac {\partial z_i^{(l+1)}}{\partial W_{i,j}^{(l)}} \\ &=\delta_i^{l+1}\cdot a_j^l\end{align} \]

原文中跳過了這一步，基於此，才能得到誤差函數對參數$\ W, b\ $每一項的偏導數。

\[\frac {\partial J_{W,b;x,y}}{\partial W_{i,j}^{(l)}} = \delta_i^{l+1}\cdot a_j^l \]

\[\frac {\partial J_{W,b;x,y}}{\partial b_{i}^{(l)}} = \delta_i^{l+1} \]

結合稀疏自編碼問題

前面的反向傳導是神經網絡的一般形式，具體到稀疏自編碼，有兩個關鍵字“稀疏”、“自”。“自”是指目標值$\ y\ \(和輸入向量\)\ x\ \(相等，“稀疏”是指過程要使隱藏層L1的大部分激活值分量接近0，個別分量顯著大於0.所以稀疏自編碼的代價函數中\)\ y\ \(直接使用\)\ x\ $的值。同時加上稀疏性懲罰項,詳見稀疏編碼一節。

\[J_sparse(W,b)=J(W,b)+\beta\sum_{j=1}^{S_2}KL(\rho||\hat\rho_j) \]

懲罰項只施加在隱藏層上，對偏導項做對應項添加，也只有隱藏層的偏導增加一項。

\[\delta_i^{(2)} = (\sum_{j=1}^{s_2}W_{ij}^{(2)}\delta_j^{(3)}+\beta (- \frac \rho {\hat \rho_i} + \frac {1-\rho} {1-\hat \rho_i}) ) f'( z_i^{(2)}) \]

后向傳導以及對網絡權值$\ W, b\ $的梯度公式不變。

在反向傳導的前因后果段落，我們已經知道了要用梯度下降法搜索使代價函數最小的自變量$\ W, b\ \(。實際上梯度下降也不是簡單減梯度，而是有下降速率參數\)\ \alpha\ $的L-BFGS，又牽扯到一個專門的優化。為了簡化要做的工作量，稀疏自編碼練習中已經幫我們集成了一個第三方實現，在minFunc文件夾中，我們只需要提供代價函數、梯度計算函數就可以調用minFunc實現梯度下降，得到最優化參數。后續的工作也就是要只需要補上sampleIMAGES.m, sparseAutoencoderCost.m, computeNumericalGradient.m的"YOUR CODE HERE"的代碼段。

exercise代碼實現難點

UFLDL給大家的學習模式很到家，把周邊的結構性代碼都寫好了matlab代碼與注釋，盡量給學習者減負。系數自編碼中主m文件是train.m。先用實現好的代價、梯度模塊調用梯度檢驗，然后將上述代價、梯度函數傳入梯度下降優化minFunc。滿足迭代次數后退出，得到訓練好的神經網絡系數矩陣，補全全部待實現模塊的完整代碼見此，https://github.com/codgeek/deeplearning.
其中主要過程是

訓練數據

訓練數據來源是圖片，第一步要在sampleIMAGES.m中將其轉化為神經網絡$\ x\ $向量。先看一下訓練數據的樣子吧

sampleIMAGES.m模塊是可以被后續課程復用的，盡量寫的精簡通用些。從一幅圖上取patchsize大小的子矩陣，要避免每個圖上取的位置相同，也要考慮並行化，循環每張圖片取patch效率較低。下文給出的方法是先隨機確定行起始行列坐標，計算出numpatches個這樣的坐標值，然后就能每個patch不關聯地取出來，才能運用parfor做並行循環。

sampleIMAGES.m
function patches = sampleIMAGES(patchsize, numpatches)
	load IMAGES;    % load images from disk dataSet=IMAGES;
    dataSet = IMAGES;
	patches = zeros(patchsize*patchsize, numpatches); % 初始化數組大小，為並行循環開辟空間
	[row, colum, numPic] = size(dataSet);
	rowStart = randi(row - patchsize + 1, numpatches, 1);% 從一幅圖上取patchsize大小子矩陣，只需要確定隨機行起始點，范圍是[1,row - patchsize + 1]
	columStart = randi(colum - patchsize + 1, numpatches, 1);% 確定隨機列起始點，范圍是colum - patchsize + 1
	randIdx = randperm(numPic); % 確定從哪一張圖上取子矩陣，打亂排列，防止生成的patch順序排列。
	parfor r=1:numpatches % 確定了起始坐標后，每個patch不關聯，可以並行循環取
    	patches(:,r) = reshape(dataSet(rowStart(r):rowStart(r) + patchsize - 1, columStart(r):columStart(r) + patchsize - 1, randIdx(floor((r-1)*numPic/numpatches)+1)),[],1);
	end
	patches = normalizeData(patches)
end

后向傳播模型的矩陣實現

稀疏自編碼最重要模塊是計算代價、梯度矩陣的sparseAutoencoderCost.m。輸入$\ W, b\ \(與訓練數據data，外加稀疏性因子、稀疏項系數、權重\)\ W\ \(系數，值得關注的是后向傳導時隱藏層輸出的殘差delta2是由隱藏層與輸出層的網絡參數\)\ W2\ \(和輸出層的殘差delta3決定的，不是輸入層與隱藏層的網絡參數\)\ W1\ $，這里最開始寫錯了，耽誤了一些時間才調試正確。下文結合代碼用注釋的形式解釋每一步具體作用。
sparseAutoencoderCost.m
function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ...
lambda, sparsityParam, beta, data)
%% 一維向量重組為便於矩陣計算的神經網絡系數矩陣。
W1 = reshape(theta(1:hiddenSizevisibleSize), hiddenSize, visibleSize);
W2 = reshape(theta(hiddenSizevisibleSize+1:2hiddenSizevisibleSize), visibleSize, hiddenSize);
b1 = theta(2hiddenSizevisibleSize+1:2hiddenSizevisibleSize+hiddenSize);
b2 = theta(2hiddenSizevisibleSize+hiddenSize+1:end);

    %% 前向傳導，W1矩陣為hiddenSize×visibleSize，訓練數據data為visibleSize×N_samples，訓練所有數據的過程正好是矩陣相乘$\ W1*data\ $。注意所有訓練數據都共用系數$\ b\ $，而單個向量的每個分量對用使用$\ b\ $的對應分量，b1*ones(1,m)是將列向量復制m遍，組成和$\ W1*data\ $相同維度的矩陣。
    [~, m] = size(data); % visibleSize×N_samples， m=N_samples
    a2 = sigmoid(W1*data + b1*ones(1,m));% active value of hiddenlayer: hiddenSize×N_samples
    a3 = sigmoid(W2*a2 + b2*ones(1,m));% output result: visibleSize×N_samples
    diff = a3 - data; % 自編碼也就意味着將激活值和原始訓練數據做差值。
    penalty = mean(a2, 2); %  measure of hiddenlayer active: hiddenSize×1
    residualPenalty =  (-sparsityParam./penalty + (1 - sparsityParam)./(1 - penalty)).*beta; % penalty factor in residual error delta2
    % size(residualPenalty)
    cost = sum(sum((diff.*diff)))./(2*m) + ...
        (sum(sum(W1.*W1)) + sum(sum(W2.*W2))).*lambda./2 + ...
        beta.*sum(KLdivergence(sparsityParam, penalty));
   
    % 后向傳導過程，隱藏層殘差需要考慮稀疏性懲罰項，公式比較清晰。
     delta3 = -(data-a3).*(a3.*(1-a3)); %  visibleSize×N_samples
     delta2 = (W2'*delta3 + residualPenalty*ones(1, m)).*(a2.*(1-a2)); %  hiddenSize×N_samples. !!! => W2'*delta3 not W1'*delta3
     % 前面已經推導出代價函數對W2的偏導，矩陣乘法里包含了公式中l層激活值a向量與1+1層殘差delta向量的點乘。
     W2grad = delta3*a2'; % ▽J(L)=delta(L+1,i)*a(l,j). sum of grade value from N_samples is got by matrix product visibleSize×N_samples * N_samples× hiddenSize. so mean value is caculated by "/N_samples"
     W1grad = delta2*data';% matrix product  visibleSize×N_samples * N_samples×hiddenSize
     
     b1grad = sum(delta2, 2);
     b2grad = sum(delta3, 2);

    % 對m個訓練數據取平均
    W1grad=W1grad./m + lambda.*W1;
    W2grad=W2grad./m + lambda.*W2;
    b1grad=b1grad./m;
    b2grad=b2grad./m;% mean value across N_sample: visibleSize ×1

    % 矩陣轉列向量
    grad = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)];

end

最終結果

在tarin.m中，將sparseAutoencoderCost.m傳入梯度優化函數minFunc，經過迭代訓練出網絡參數，通常見到的各個方向的邊緣檢測圖表示的是權重矩陣對每個隱藏層激活值的特征，當輸入值為對應特征值時，每個激活值會有最大響應，同時的其余隱藏節點處於抑制狀態。所以只需要把W1矩陣每一行的64個向量還原成8*8的圖片patch，也就是特征值了，每個隱藏層對應一個，總共100個。結果如下圖。

接下來我們驗證一下隱藏層對應的輸出是否滿足激活、抑制的要求，將上述輸入值用參數矩陣傳導到隱藏層，也就是$\ W1*W1'\ $.可見，每個輸入對應的隱藏層只有一個像素是白色，表示接近1，其余是暗色調接近0，滿足了稀疏性要求，而且激活的像素位置是順序排列的。

繼續改變不同的輸入層單元個數，隱藏層單元個數，可以得到更多有意思的結果！