前言

實驗環境：win7， matlab2015b，16G記憶體，2T機械硬碟

難點：本實驗難點在於執行時間比較長，跑一次都快一天了，並且我還要驗證各種代價函式的對錯，所以跑了很多次。

實驗基礎說明：

     ①不同點：本節實驗中的基是標準正交的，也是線性獨立的，而Deep Learning 12_深度學習UFLDL教程：Sparse Coding_exercise（斯坦福大學深度學習教程）中的基的個數比原始資料還多，所以它的基不一定是線性獨立的，儘管稀疏編碼學習到的基在某些情況下這已經滿足需要，但有時我們仍然希望得到的是一組線性獨立基。在稀疏編碼中，權值矩陣 A 用於將特徵 s

     ②相同點：稀疏編碼和ICA演算法學習到的特徵都要求是稀疏的。

2.Ng為什麼說投影步驟實際上可以看成就是ZCA白化？下面是我的個人理解

通過主成分分析和白化可知，ZCA白化整個過程實際上是：

而投影步驟的整個過程是：

通過對比以上兩者公式可知，投影步驟可以看成是ε為0的ZCA白化（即：無正則ZCA白化），其中，ZCA白化中的資料旋轉方向U相當於投影步驟中的W，ZCA白化中的特徵值λ相當於投影步驟中的(WW^T

3.資料必須經過無正則 ZCA白化（也即,ε設為0），但是為什麼？下面是我的個人理解

因為前面已經說了，標準正交基W的投影步驟實際上就是無正則ZCA白化，而特徵為Wx，用特徵和基表示的原始資料為W^TWx（原因見下面的解釋），為了統一，使特徵Wx和W^TWx也是白化的，就需要也對原始資料x做無正則ZCA白化

4.代價函式及其梯度

標準約束項WW^T = I通過投影步驟實現，所以實現時代價函式為：

其代價函式對W的梯度或偏導為：      

5.代價函式的推導

      通過以前本人對Ng的代價函式推導的理解（即：Deep Learning 12_深度學習UFLDL教程：Sparse Coding_exercise（斯坦福大學深度學習教程）），推導本實驗的代價函式，而Deep learning：三十九(ICA模型練習)中並沒有推導這些，而且從我的推導結果可知他的代價函式形式是錯誤的（準確地說，可能也不叫錯誤，因為它只是增加了一項非常複雜的多餘的項，這一項的值永遠不會變，並且使執行時間加大了1倍多，但還是可以提取出標準正交基），我覺得這是因為他並沒有理解代價函式是怎麼來的，實際這一點完全可從他整篇文章中看出來。

 本節實驗中，我們希望學習得到一組基向量――以列向量形式構成的矩陣 W，其滿足以下特點：首先，與稀疏編碼一樣，特徵是稀疏的；其次，基是標準正交的。

①稀疏懲罰項：因為要求學習到的特徵是稀疏的，且學習到的特徵表示為Wx，所以代價函式必須要有如下稀疏懲罰項：

                                                ，                          其中，是x的Lk範數，等價於 。L1 範數是向量元素的絕對值之和

又因為以後我們為求代價函式的最小值，會對代價函式求導，而對s在0點處不可導的（理解：y=|x|在x=0處不可導），所以為了方便以後求代價函式最小值，可把替換為，其中ε 是“平滑引數”（"smoothing parameter"）或者“稀疏引數”（"sparsity parameter"） 。所以把以上代價函式改為：

其中，m為樣本個數，λ為1/m，是稀疏懲罰項的權重。

②標準正交項：因為要求學習到的基W是標準正交的，所以代價函式必須有以下約束：

                                                                                                              WW^T = I

      在實際實現這一約束過程中，並沒有把代價函式中應包含的這一標準正交約束加入到代價函式中，而是把它放到利用梯度下降法優化代價函式階段。也就是，用梯度下降法優化代價函式時，在梯度下降的每一步中增加投影步驟，以滿足標準正交約束。所以，在程式設計寫代價函式並求它梯度（即：orthonormalICACost.m）時，實際上並沒有包含這一標準正交項。

③重構項：因為本節實驗的目的是尋找資料X的標準正交基W，並把它顯示出來。而學習到的特徵表示為Wx，所以聯絡上節實驗（即：Deep Learning 12_深度學習UFLDL教程：Sparse Coding_exercise（斯坦福大學深度學習教程））可知把資料x用標準正交基W表示出來的形式為：W^TWx，即：x=基*特徵集 。

     為了使資料x和W^TWx之間的誤差最小，那麼代價函式必須要包括這兩者的均方差，並且要使這個均方差最小，即最小化如下項：

                                                                                                                                  ^{，其中，L2範數就是平方和的開方。}

上面的形式只是Adrew Ng在用反向傳導思想求導中的表達，實際上需要再除以樣本數m才是真正的均方差。

所以代價函式也應該包含上面這個重構項。但是因為代價函式有標準正交約束WW^T = I，所以可推知只要W滿足標準正交約束，那麼這個重構項就永遠為0。這一點並不僅是我猜測，我還通過程式設計實踐證明這個重構項恆等於0.003550。證明很簡單，只需要把Deep learning：三十九(ICA模型練習)中orthonormalICACost.m改為如下即可：

function [cost, grad] = orthonormalICACost(theta, visibleSize, numFeatures, patches, epsilon)
%orthonormalICACost - compute the cost and gradients for orthonormal ICA
%                     (i.e. compute the cost ||Wx||_1 and its gradient)

    weightMatrix = reshape(theta, numFeatures, visibleSize);
    
    cost = 0;
    grad = zeros(numFeatures, visibleSize);
    
    % -------------------- YOUR CODE HERE --------------------
    % Instructions:
    %   Write code to compute the cost and gradient with respect to the
    %   weights given in weightMatrix.     
    % -------------------- YOUR CODE HERE --------------------     

%% 方法
      lambda =  8e-6;% 0.5e-4
    num_samples = size(patches,2);
    
    cost_part1 = sum(sum((weightMatrix'*weightMatrix*patches-patches).^2))./num_samples;
    cost_part2 = sum(sum(sqrt((weightMatrix*patches).^2+epsilon)))*lambda;
    cost = cost_part1  +  cost_part2;
    grad = (2*weightMatrix*(weightMatrix'*weightMatrix*patches-patches)*patches'+...
        2*weightMatrix*patches*(weightMatrix'*weightMatrix*patches-patches)')./num_samples+...
        (weightMatrix*patches./sqrt((weightMatrix*patches).^2+epsilon))*patches'*lambda;
    
    grad = grad(:);
    fprintf('%11s%16s\n','cost_part1','cost_part2');
    fprintf('   %14.6f  %14.6f\n', cost_part1, cost_part2);
    
end

然後執行ICAExercise.m，執行結果中cost_part1就是重構項的值，它恆等於0.003550，不會隨著迭代次數的增加而減少。

因此在orthonormalICACost.m中實現的代價函式，實際上只包含了稀疏懲罰項。

綜上，orthonormalICACost.m中代價函式為：

實驗步驟

1.初始化引數，載入資料庫。

2.把資料進行無正則ZCA白化。注意把ε設為0。

3.實現ICA演算法的代價函式及梯度計算（見orthonormalICACost.m），並檢查梯度計算是否正確。注意，這一步中並沒有把代價函式中應包含的標準正交約束加入到代價函式中，而是把它放下一步，即利用梯度下降法優化代價函式階段。具體做法就是：在梯度下降的每一步中增加投影步驟，以滿足標準正交約束。

4.利用梯度下降法優化代價函式，同時在梯度下降的每一步中增加投影步驟，以滿足標準正交約束，並且用線搜尋演算法（見Convex Optimization by Boyd and Vandenbergh）來加速梯度。

5.得到優化結果，並把權值矩陣weightMatrix（它就是我們要求的線性獨立標準正交基）顯示出來。

結果

 部分原始資料：

迭代1萬次的結果：

 迭代2萬次的結果：

 迭代5萬次的結果：

對比可知，5萬次的執行結果比Ng的結果更好一點。

我也運行了Deep learning：三十九(ICA模型練習)中的程式碼，實現了作者的意圖，但是他的程式碼實在執行時間太長。整個實驗實際上不算難，主要是程式碼執行時間長，以及代價函式的理解。

程式碼

 ICAExercise.m

%% CS294A/CS294W Independent Component Analysis (ICA) Exercise

%  Instructions
%  ------------
% 
%  This file contains code that helps you get started on the
%  ICA exercise. In this exercise, you will need to modify
%  orthonormalICACost.m and a small part of this file, ICAExercise.m.

%%======================================================================
%% STEP 0: Initialization
%  Here we initialize some parameters used for the exercise.

numPatches = 20000;
numFeatures = 121;
imageChannels = 3;
patchDim = 8;
visibleSize = patchDim * patchDim * imageChannels;

outputDir = '.';
epsilon = 1e-6; % L1-regularisation epsilon |Wx| ~ sqrt((Wx).^2 + epsilon)

%%======================================================================
%% STEP 1: Sample patches

patches = load('stlSampledPatches.mat');
patches = patches.patches(:, 1:numPatches);
displayColorNetwork(patches(:, 1:100));

%%======================================================================
%% STEP 2: ZCA whiten patches
%  In this step, we ZCA whiten the sampled patches. This is necessary for
%  orthonormal ICA to work.

patches = patches / 255;
meanPatch = mean(patches, 2);
patches = bsxfun(@minus, patches, meanPatch);

sigma = patches * patches';
[u, s, v] = svd(sigma);
ZCAWhite = u * diag(1 ./ sqrt(diag(s))) * u';
patches = ZCAWhite * patches;

%%======================================================================
%% STEP 3: ICA cost functions
%  Implement the cost function for orthornomal ICA (you don't have to 
%  enforce the orthonormality constraint in the cost function) 
%  in the function orthonormalICACost in orthonormalICACost.m.
%  Once you have implemented the function, check the gradient.

% Use less features and smaller patches for speed
debug = false;
if debug
numFeatures = 5;
patches = patches(1:3, 1:5);
visibleSize = 3;
numPatches = 5;

weightMatrix = rand(numFeatures, visibleSize);

[cost, grad] = orthonormalICACost(weightMatrix, visibleSize, numFeatures, patches, epsilon);

numGrad = computeNumericalGradient( @(x) orthonormalICACost(x, visibleSize, numFeatures, patches, epsilon), weightMatrix(:) );
% Uncomment to display the numeric and analytic gradients side-by-side
% disp([numGrad grad]); 
diff = norm(numGrad-grad)/norm(numGrad+grad);
fprintf('Orthonormal ICA difference: %g\n', diff);
assert(diff < 1e-7, 'Difference too large. Check your analytic gradients.');

fprintf('Congratulations! Your gradients seem okay.\n');
end
%%======================================================================
%% STEP 4: Optimization for orthonormal ICA
%  Optimize for the orthonormal ICA objective, enforcing the orthonormality
%  constraint. Code has been provided to do the gradient descent with a
%  backtracking line search using the orthonormalICACost function 
%  (for more information about backtracking line search, you can read the 
%  appendix of the exercise).
%
%  However, you will need to write code to enforce the orthonormality 
%  constraint by projecting weightMatrix back into the space of matrices 
%  satisfying WW^T  = I.
%
%  Once you are done, you can run the code. 10000 iterations of gradient
%  descent will take around 2 hours, and only a few bases will be
%  completely learned within 10000 iterations. This highlights one of the
%  weaknesses of orthonormal ICA - it is difficult to optimize for the
%  objective function while enforcing the orthonormality constraint - 
%  convergence using gradient descent and projection is very slow.

weightMatrix = rand(numFeatures, visibleSize);

[cost, grad] = orthonormalICACost(weightMatrix(:), visibleSize, numFeatures, patches, epsilon);

fprintf('%11s%16s%10s\n','Iteration','Cost','t');

startTime = tic();

% Initialize some parameters for the backtracking line search
alpha = 0.5;
t = 0.02;
lastCost = 1e40;

% Do 10000 iterations of gradient descent
for iteration = 1:50000
                       
    grad = reshape(grad, size(weightMatrix));
    newCost = Inf;        
    linearDelta = sum(sum(grad .* grad));
    
    % Perform the backtracking line search
    while 1
        considerWeightMatrix = weightMatrix - alpha * grad;
        % -------------------- YOUR CODE HERE --------------------
        % Instructions:
        %   Write code to project considerWeightMatrix back into the space
        %   of matrices satisfying WW^T = I.
        %   
        %   Once that is done, verify that your projection is correct by 
        %   using the checking code below. After you have verified your
        %   code, comment out the checking code before running the
        %   optimization.
        
%         % Project considerWeightMatrix such that it satisfies WW^T = I
%         error('Fill in the code for the projection here');        
        considerWeightMatrix = (considerWeightMatrix*considerWeightMatrix')^(-0.5)*considerWeightMatrix;
        % Verify that the projection is correct
        temp = considerWeightMatrix * considerWeightMatrix';
        temp = temp - eye(numFeatures);
        assert(sum(temp(:).^2) < 1e-23, 'considerWeightMatrix does not satisfy WW^T = I. Check your projection again');
        error('Projection seems okay. Comment out verification code before running optimization.');
        
        % -------------------- YOUR CODE HERE --------------------                                        

        [newCost, newGrad] = orthonormalICACost(considerWeightMatrix(:), visibleSize, numFeatures, patches, epsilon);
        if newCost > lastCost - alpha * t * linearDelta
%             fprintf('   %14.6f  %14.6f\n', newCost, lastCost - alpha * t * linearDelta);
            t = 0.9 * t;
        else
            break;
        end
    end
   
    lastCost = newCost;
    weightMatrix = considerWeightMatrix;
    
    fprintf('  %9d  %14.6f  %8.7g\n', iteration, newCost, t);
    
    t = 1.1 * t;
    
    cost = newCost;
    grad = newGrad;
           
    % Visualize the learned bases as we go along    
    if mod(iteration, 1000) == 0
        duration = toc(startTime);
        % Visualize the learned bases over time in different figures so 
        % we can get a feel for the slow rate of convergence
        figure(floor(iteration /  1000));
        displayColorNetwork(weightMatrix'); 
    end
                   
end

% Visualize the learned bases
displayColorNetwork(weightMatrix');

orthonormalICACost.m

function [cost, grad] = orthonormalICACost(theta, visibleSize, numFeatures, patches, epsilon)
%orthonormalICACost - compute the cost and gradients for orthonormal ICA
%                     (i.e. compute the cost ||Wx||_1 and its gradient)

    weightMatrix = reshape(theta, numFeatures, visibleSize);
    
    cost = 0;
    grad = zeros(numFeatures, visibleSize);
    
    % -------------------- YOUR CODE HERE --------------------
    % Instructions:
    %   Write code to compute the cost and gradient with respect to the
    %   weights given in weightMatrix.     
    % -------------------- YOUR CODE HERE --------------------     

%% 
    num_samples = size(patches,2); %樣本個數

    aux1 = sqrt(((weightMatrix*patches).^2) + epsilon);
    cost = sum(aux1(:))/num_samples;
    grad = ((weightMatrix*patches)./aux1)*patches'./num_samples;
    grad = grad(:);

    
end

參考資料