稀疏自編碼器的學習結構:
稀疏自編碼器Ⅰ：
神經網絡
反向傳導算法
梯度檢驗與高級優化
稀疏自編碼器Ⅱ：
自編碼算法與稀疏性
可視化自編碼器訓練結果
Exercise: Sparse Autoencoder
自編碼算法與稀疏性
已經討論了神經網絡在有監督學習中的應用，其中訓練樣本是有類別標簽的（x_i，y_i）。
自編碼神經網絡是一種無監督學習算法，它使用了反向傳播算法，並讓目標值等於輸入值x_i = y_i 。
下圖是一個自編碼神經網絡的示例。
一次autoencoder學習，結構三層：輸入層單元數=輸出層單元數，隱藏層。自編碼神經網絡嘗試學習一個輸入約等於輸出的恆等函數。
給自編碼神經網絡加入某些限制，我們就可以從輸入數據中發現一些有趣的結構。比如限定隱藏神經元的數量。
隱藏層神經元的數量較小
比如使隱藏層神經元的數量<輸入層神經元的數量<輸出層經元的數量，迫使自編碼神經網絡去學習輸入數據的壓縮表示。
當輸入數據是完全隨機的，比如輸入特征完全無關，難學習。
當輸入數據中隱含着一些特定的結構，比如輸入特征是彼此相關的，算法就可以發現輸入數據中的這些相關性。
事實上，這一簡單的自編碼神經網絡通常可以學習出一個跟主元分析（PCA）結果非常相似的輸入數據的低維表示。
隱藏層神經元的數量較大
仍然通過給自編碼神經網絡施加一些其他的限制條件來發現輸入數據中的結構。
給隱藏層神經元加入稀疏性限制
稀疏性可以被簡單地解釋如下。如果當神經元的輸出接近於1的時候我們認為它被激活，而輸出接近於0的時候認為它被抑制，那么使得神經元大部分的時間都是被抑制的限制則被稱作稀疏性限制。這里我們假設的神經元的激活函數是sigmoid函數。如果你使用tanh作為激活函數的話，當神經元輸出為-1的時候，我們認為神經元是被抑制的。
以上是稀疏性的含義，具體獲得稀疏性的方法ufldl教程中有詳細講述，這里只說核心概念框架。
隱藏神經元 j的平均活躍度（在訓練集上取平均）
注意，計算用到了前向傳播算法，而BP也用到了，內存夠保存算一遍，否則兩遍。
限制
其中p是稀疏性參數，通常是一個接近於0的較小的值（比如 p=0.05 ）
為了實現這一限制，我們將會在我們的優化目標函數中加入一個額外的懲罰因子，而這一懲罰因子將懲罰那些 和 有顯著不同的情況從而使得隱藏神經元的平均活躍度保持在較小范圍內（稀疏性）。
懲罰因子的具體形式有很多種合理的選擇，我們將會選擇以下這一種：
KL divergence 性質：相等為0，隨着之間的差異增大而單調遞增。
現在，我們的總體代價函數可以表示為：
To incorporate the KL-divergence term into your derivative calculation, there is a simple-to-implement trick involving only a small change to your code.
具體地，在BP第二層中：
可視化自編碼器訓練結果
訓練完（稀疏）自編碼器，我們還想把這自編碼器學到的函數可視化出來，好弄明白它到底學到了什么。我們以在10×10圖像（即n=100）上訓練自編碼器為例。
Exercise: Sparse Autoencoder實驗
首先是從如下圖這樣的10 image（512×512），中sample 10000 image patches（8×8）。
sample 10000 image patches（8×8） and concatenate them into a 64×10000 matrix
display a random sample of 204 patches from the dataset
以下是一些設置參數，就是學習單個自動編碼器，學到一個隱含層的基（參數W，b）。
| visibleSize = 8*8; % number of input units hiddenSize = 25; % number of hidden units sparsityParam = 0.01; % 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 = 0.0001; % weight decay parameter beta = 3; % weight of sparsity penalty term |
以下是序列輸出結果
| >> train Iteration FunEvals Step Length Function Val Opt Cond 1 3 8.63782e-03 3.98056e+01 1.03759e+03 2 4 1.00000e+00 6.68382e+00 2.49253e+02 ……. 399 414 1.00000e+00 4.49948e-01 1.45238e-02 400 415 1.00000e+00 4.49947e-01 1.40765e-02 Exceeded Maximum Number of Iterations 時間已過 20.942873 秒。 |
教程中說Our implementation took around 5 minutes to run on a fast computer.
sparse autoencoder algorithm learning a set of edge detectors.
隱藏層權重可視化后，我們可以看出學到了一組邊緣檢測器，一組基或稱字典。
練習代碼：
train
%% CS294A/CS294W Programming Assignment Starter Code % Instructions % ------------ % % This file contains code that helps you get started on the % programming assignment. You will need to complete the code in sampleIMAGES.m, % sparseAutoencoderCost.m and computeNumericalGradient.m. % For the purpose of completing the assignment, you do not need to % change the code in this file. % %%====================================================================== %% STEP 0: Here we provide the relevant parameters values that will % allow your sparse autoencoder to get good filters; you do not need to % change the parameters below. visibleSize = 8*8; % number of input units hiddenSize = 25; % number of hidden units sparsityParam = 0.01; % 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 = 0.0001; % weight decay parameter beta = 3; % weight of sparsity penalty term %%====================================================================== %% STEP 1: Implement sampleIMAGES % % After implementing sampleIMAGES, the display_network command should % display a random sample of 200 patches from the dataset patches = sampleIMAGES; display_network(patches(:,randi(size(patches,2),204,1)),8); % Obtain random parameters theta theta = initializeParameters(hiddenSize, visibleSize); %%====================================================================== %% STEP 2: Implement sparseAutoencoderCost % % You can implement all of the components (squared error cost, weight decay term, % sparsity penalty) in the cost function at once, but it may be easier to do % it step-by-step and run gradient checking (see STEP 3) after each step. We % suggest implementing the sparseAutoencoderCost function using the following steps: % % (a) Implement forward propagation in your neural network, and implement the % squared error term of the cost function. Implement backpropagation to % compute the derivatives. Then (using lambda=beta=0), run Gradient Checking % to verify that the calculations corresponding to the squared error cost % term are correct. % % (b) Add in the weight decay term (in both the cost function and the derivative % calculations), then re-run Gradient Checking to verify correctness. % % (c) Add in the sparsity penalty term, then re-run Gradient Checking to % verify correctness. % % Feel free to change the training settings when debugging your % code. (For example, reducing the training set size or % number of hidden units may make your code run faster; and setting beta % and/or lambda to zero may be helpful for debugging.) However, in your % final submission of the visualized weights, please use parameters we % gave in Step 0 above. [cost, grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, lambda, ... sparsityParam, beta, patches); %%====================================================================== %% STEP 3: Gradient Checking % % Hint: If you are debugging your code, performing gradient checking on smaller models % and smaller training sets (e.g., using only 10 training examples and 1-2 hidden % units) may speed things up. % First, lets make sure your numerical gradient computation is correct for a % simple function. After you have implemented computeNumericalGradient.m, % run the following: % checkNumericalGradient(); % % % Now we can use it to check your cost function and derivative calculations % % for the sparse autoencoder. % numgrad = computeNumericalGradient( @(x) sparseAutoencoderCost(x, visibleSize, ... % hiddenSize, lambda, ... % sparsityParam, beta, ... % patches), theta); % % % Use this to visually compare the gradients side by side % disp([numgrad grad]); % % % Compare numerically computed gradients with the ones obtained from backpropagation % diff = norm(numgrad-grad)/norm(numgrad+grad); % disp(diff); % Should be small. In our implementation, these values are % % usually less than 1e-9. % % % When you got this working, Congratulations!!! %%====================================================================== %% STEP 4: After verifying that your implementation of % sparseAutoencoderCost is correct, You can start training your sparse % autoencoder with minFunc (L-BFGS). % Randomly initialize the parameters theta = initializeParameters(hiddenSize, visibleSize); % Use minFunc to minimize the function 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'; [opttheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ... visibleSize, hiddenSize, ... lambda, sparsityParam, ... beta, patches), ... theta, options); %%====================================================================== %% STEP 5: Visualization W1 = reshape(opttheta(1:hiddenSize*visibleSize), hiddenSize, visibleSize); display_network(W1', 12); print -djpeg weights.jpg % save the visualization to a file
sampleIMAGES
function patches = sampleIMAGES() % sampleIMAGES % Returns 10000 patches for training load IMAGES; % load images from disk patchsize = 8; % we'll use 8x8 patches numpatches = 10000; % Initialize patches with zeros. Your code will fill in this matrix--one % column per patch, 10000 columns. patches = zeros(patchsize*patchsize, numpatches); %% ---------- YOUR CODE HERE -------------------------------------- % Instructions: Fill in the variable called "patches" using data % from IMAGES. % % IMAGES is a 3D array containing 10 images % For instance, IMAGES(:,:,6) is a 512x512 array containing the 6th image, % and you can type "imagesc(IMAGES(:,:,6)), colormap gray;" to visualize % it. (The contrast on these images look a bit off because they have % been preprocessed using using "whitening." See the lecture notes for % more details.) As a second example, IMAGES(21:30,21:30,1) is an image % patch corresponding to the pixels in the block (21,21) to (30,30) of % Image 1 for imageNum = 1 : 10 image = IMAGES(:, :, imageNum); [rowNum, colNum] = size(image); for patchNum = 1 : 1000 xPos = randi(rowNum - patchsize + 1); yPos = randi(colNum - patchsize + 1); patches(:, 1000 * (imageNum - 1) + patchNum) = ... reshape(image(xPos : xPos + 7, yPos : yPos + 7), 64, 1); end end %% --------------------------------------------------------------- % For the autoencoder to work well we need to normalize the data % Specifically, since the output of the network is bounded between [0,1] % (due to the sigmoid activation function), we have to make sure % the range of pixel values is also bounded between [0,1] patches = normalizeData(patches); end %% --------------------------------------------------------------- function patches = normalizeData(patches) % Squash data to [0.1, 0.9] since we use sigmoid as the activation % function in the output layer % Remove DC (mean of images). patches = bsxfun(@minus, patches, mean(patches)); % Truncate to +/-3 standard deviations and scale to -1 to 1 pstd = 3 * std(patches(:)); patches = max(min(patches, pstd), -pstd) / pstd; % Rescale from [-1,1] to [0.1,0.9] patches = (patches + 1) * 0.4 + 0.1; end
sparseAutoencoderCost
function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ... lambda, sparsityParam, beta, data) % visibleSize: the number of input units (probably 64) % hiddenSize: the number of hidden units (probably 25) % lambda: weight decay parameter % sparsityParam: The desired average activation for the hidden units (denoted in the lecture % notes by the greek alphabet rho, which looks like a lower-case "p"). % beta: weight of sparsity penalty term % data: Our 64x10000 matrix containing the training data. So, data(:,i) is the i-th training example. % The input theta is a vector (because minFunc expects the parameters to be a vector). % We first convert theta to the (W1, W2, b1, b2) matrix/vector format, so that this % follows the notation convention of the lecture notes. W1 = reshape(theta(1:hiddenSize*visibleSize), hiddenSize, visibleSize); W2 = reshape(theta(hiddenSize*visibleSize+1:2*hiddenSize*visibleSize), visibleSize, hiddenSize); b1 = theta(2*hiddenSize*visibleSize+1:2*hiddenSize*visibleSize+hiddenSize); b2 = theta(2*hiddenSize*visibleSize+hiddenSize+1:end); % Cost and gradient variables (your code needs to compute these values). % Here, we initialize them to zeros. cost = 0; W1grad = zeros(size(W1)); W2grad = zeros(size(W2)); b1grad = zeros(size(b1)); b2grad = zeros(size(b2)); %% ---------- YOUR CODE HERE -------------------------------------- % 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. % Jcost = 0; Jweight = 0; Jsparse = 0; [n m] = size(data); z2 = W1*data+repmat(b1,1,m); a2 = sigmoid(z2); z3 = W2*a2+repmat(b2,1,m); a3 = sigmoid(z3); Jcost = (0.5/m)*sum(sum((a3-data).^2)); Jweight = (1/2)*(sum(sum(W1.^2))+sum(sum(W2.^2))); rho = (1/m).*sum(a2,2); Jsparse = sum(sparsityParam.*log(sparsityParam./rho)+ ... (1-sparsityParam).*log((1-sparsityParam)./(1-rho))); cost = Jcost+lambda*Jweight+beta*Jsparse; d3 = -(data-a3).*dsigmoid(a3); sterm = beta*(-sparsityParam./rho+(1-sparsityParam)./(1-rho)); d2 = (W2'*d3+repmat(sterm,1,m)).*dsigmoid(a2); W1grad = W1grad+d2*data'; W1grad = (1/m)*W1grad+lambda*W1; W2grad = W2grad+d3*a2'; W2grad = (1/m).*W2grad+lambda*W2; b1grad = b1grad+sum(d2,2); b1grad = (1/m)*b1grad; b2grad = b2grad+sum(d3,2); b2grad = (1/m)*b2grad; %------------------------------------------------------------------- % 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 = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)]; 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)). function sigm = sigmoid(x) sigm = 1 ./ (1 + exp(-x)); end function dsigm = dsigmoid(a) dsigm = a .* (1.0 - a); end
computeNumericalGradient
function numgrad = computeNumericalGradient(J, theta) % numgrad = computeNumericalGradient(J, theta) % theta: a vector of parameters % J: a function that outputs a real-number. Calling y = J(theta) will return the % function value at theta. % Initialize numgrad with zeros numgrad = zeros(size(theta)); %% ---------- YOUR CODE HERE -------------------------------------- % Instructions: % Implement numerical gradient checking, and return the result in numgrad. % (See Section 2.3 of the lecture notes.) % You should write code so that numgrad(i) is (the numerical approximation to) the % partial derivative of J with respect to the i-th input argument, evaluated at theta. % I.e., numgrad(i) should be the (approximately) the partial derivative of J with % respect to theta(i). % % Hint: You will probably want to compute the elements of numgrad one at a time. epsilon = 1e-4; n = size(theta,1); E = eye(n); for i = 1:n delta = E(:,i)*epsilon; numgrad(i) = (J(theta+delta)-J(theta-delta))/(epsilon*2.0); end %% --------------------------------------------------------------- end
聲明：
不要將本博客用作為商業目的，並且本人博客使用的內容僅僅作為個人學習，其中包括引用的或沒找到出處而未引用的內容。
轉載請注明，本文地址：http://www.cnblogs.com/JayZen/p/4129386.html
參考：
UFLDL教程
Hugo Larochelle nn course