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 = ;

    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-;% 0.5e-4

    num_samples = size(patches,);

    cost_part1 = sum(sum((weightMatrix'*weightMatrix*patches-patches).^2))./num_samples;

    cost_part2 = sum(sum(sqrt((weightMatrix*patches).^+epsilon)))*lambda;

    cost = cost_part1  +  cost_part2;

    grad = (*weightMatrix*(weightMatrix'*weightMatrix*patches-patches)*patches'+...

        *weightMatrix*patches*(weightMatrix'*weightMatrix*patches-patches)')./num_samples+...

        (weightMatrix*patches./sqrt((weightMatrix*patches).^+epsilon))*patches'*lambda;

    grad = grad(:);

    fprintf('%11s%16s\n','cost_part1','cost_part2');

    fprintf('   %14.6f  %14.6f\n', cost_part1, cost_part2);

end

%% 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 : Initialization

%  Here we initialize some parameters used for the exercise.

numPatches = ;

numFeatures = ;

imageChannels = ;

patchDim = ;

visibleSize = patchDim * patchDim * imageChannels;

outputDir = '.';

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

%%======================================================================

%% STEP : Sample patches

patches = load('stlSampledPatches.mat');

patches = patches.patches(:, :numPatches);

displayColorNetwork(patches(:, :));

%%======================================================================

%% STEP : ZCA whiten patches

%  In this step, we ZCA whiten the sampled patches. This is necessary for

%  orthonormal ICA to work.

patches = patches / ;

meanPatch = mean(patches, );

patches = bsxfun(@minus, patches, meanPatch);

sigma = patches * patches';

[u, s, v] = svd(sigma);

ZCAWhite = u * diag( ./ sqrt(diag(s))) * u';

patches = ZCAWhite * patches;

%%======================================================================

%% STEP : 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 = ;

patches = patches(:, :);

visibleSize = ;

numPatches = ;

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-, 'Difference too large. Check your analytic gradients.');

fprintf('Congratulations! Your gradients seem okay.\n');

end

%%======================================================================

%% STEP : 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.  iterations of gradient

%  descent will take around  hours, and only a few bases will be

%  completely learned within  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  iterations of gradient descent

for iteration = :

    grad = reshape(grad, size(weightMatrix));

    newCost = Inf;

    linearDelta = sum(sum(grad .* grad));

    % Perform the backtracking line search

    while

        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(:).^) < 1e-, '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, ) ==

        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 /  ));

        displayColorNetwork(weightMatrix');

    end

end

% Visualize the learned bases

displayColorNetwork(weightMatrix');

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 = ;

    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,); %样本个数

    aux1 = sqrt(((weightMatrix*patches).^) + epsilon);

    cost = sum(aux1(:))/num_samples;

    grad = ((weightMatrix*patches)./aux1)*patches'./num_samples;

    grad = grad(:);

end