稀疏编码(SparseCoding)
sparse coding也是deep learning中一个重要的分支，同样能够提取出数据集很好的特征(稀疏的)。选择使用具有稀疏性的分量来表示我们的输入数据是有原因的，因为绝大多数的感官数据，比如自然图像，可以被表示成少量基本元素的叠加，在图像中这些基本元素可以是面或者线。
稀疏编码算法的目的就是找到一组基向量
这里构成的基向量要求是超完备的，即要求k大于n，这样的方程就大多情况会有无穷多个解，此时我们给它加入一个稀疏性的限制，最终的优化公式变成了如下形式：
其中S(ai)就是稀疏惩罚项，可以是L0或者L1范数，L1范数和L0范数都可以表示稀疏编码，L1范数是L0范数的最优凸近似，但是L1因具有比L0更好的优化求解特性而被广泛应用。
稀疏编码自编码表达：
将稀疏编码用到深度学习中，用于提取数据集良好的稀疏特征，设A为超完备的基向量，s表示输入数据最后的稀疏特征(也就是稀疏编码中的稀疏系数),这样就可以表示成X = A*s。
其实这里的A就等同于稀疏自编码中的W2，而s就是隐层的结点值。(当具有很多样本的时候，s就是一个矩阵，每一列表示的是一个样本的稀疏特征值)
最终的优化函数变成了：
优化步骤为：
可以采用以下两个trick来提高最后的迭代速度和精度。
(1)将样本表示成良好的“迷你块”，比如你有10000个样本，我们可以每次只随机选择2000个mini_patches进行迭代，这样不仅提高了迭代的速度，也提升了收敛的速度。
(2)良好的s初始化。因为我们的目标是X=As,因此交叉迭代优化的过程中，在A给定的情况下，我们可以将s初始化为s=AT*X,但这样可能会导致稀疏性的缺失，我们在做一个规范，这里其中的s的行数就等于A列数，然后用s的每个元素除以其所在A的那一列向量的2范数。
即 
因此最后的优化算法步骤表示为：
注意：以上对s和A采用交叉迭代优化，其中我们会发现分别对s和A求导的时候，发现可以直接得出A的解析形式的解，因此在优化A的时候直接给出其解析形式的解就可以了，而s我们无法给出其解析形式的解，就需要用梯度迭代等无约束的优化问题了。
测试：当有一个新的样本x，我们需要用以上训练集得到的A来优化以上cost函数得到s就是该样本的稀疏特征。这样相比之前的前馈网络的，每次对新的数据样本进行“编码”，我们必须再次执行优化过程来得到所需的系数。
注：教程中拓扑稀疏编码的内容我还没有弄明白是如何得到groupMatrix的，这里就不误导大家了。
练习答案：
以下对grad的求解可以参看这篇blog：http://www.cnblogs.com/tornadomeet/archive/2013/04/16/3024292.html给出的推导结果。
Sparse_coding_exercise.m
%% CS294A/CS294W Sparse Coding Exercise
% Instructions
% ------------
%
% This file contains code that helps you get started on the
% sparse coding exercise. In this exercise, you will need to modify
% sparseCodingFeatureCost.m and sparseCodingWeightCost.m. You will also
% need to modify this file, sparseCodingExercise.m slightly.
% Add the paths to your earlier exercises if necessary
% addpath /path/to/solution
%%======================================================================
%% STEP 0: Initialization
% Here we initialize some parameters used for the exercise.
numPatches = 20000; % number of patches
numFeatures = 121; % number of features to learn
patchDim = 8; % patch dimension
visibleSize = patchDim * patchDim;
% dimension of the grouping region (poolDim x poolDim) for topographic sparse coding
poolDim = 3;
% number of patches per batch
batchNumPatches = 2000;
lambda = 5e-5; % L1-regularisation parameter (on features)
epsilon = 1e-5; % L1-regularisation epsilon |x| ~ sqrt(x^2 + epsilon)
gamma = 1e-2; % L2-regularisation parameter (on basis)
%%======================================================================
%% STEP 1: Sample patches
images = load('IMAGES.mat');
images = images.IMAGES;
patches = sampleIMAGES(images, patchDim, numPatches);
display_network(patches(:, 1:64));
%%======================================================================
%% STEP 2: Implement and check sparse coding cost functions
% Implement the two sparse coding cost functions and check your gradients.
% The two cost functions are
% 1) sparseCodingFeatureCost (in sparseCodingFeatureCost.m) for the features
% (used when optimizing for s, which is called featureMatrix in this exercise)
% 2) sparseCodingWeightCost (in sparseCodingWeightCost.m) for the weights
% (used when optimizing for A, which is called weightMatrix in this exericse)
% We reduce the number of features and number of patches for debugging
% numFeatures = 25;
% patches = patches(:, 1:5);
% numPatches = 5;
weightMatrix = randn(visibleSize, numFeatures) * 0.005;
featureMatrix = randn(numFeatures, numPatches) * 0.005;
%% STEP 2a: Implement and test weight cost
% Implement sparseCodingWeightCost in sparseCodingWeightCost.m and check
% the gradient
[cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon);
numgrad = computeNumericalGradient( @(x) sparseCodingWeightCost(x, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon), weightMatrix(:) );
% Uncomment the blow line to display the numerical and analytic gradients side by side
% disp([numgrad grad]);
diff = norm(numgrad-grad)/norm(numgrad+grad);
fprintf('Weight difference: %g\n', diff);
assert(diff < 1e-8, 'Weight difference too large. Check your weight cost function. ');
%% STEP 2b: Implement and test feature cost (non-topographic)
% Implement sparseCodingFeatureCost in sparseCodingFeatureCost.m and check
% the gradient. You may wish to implement the non-topographic version of
% the cost function first, and extend it to the topographic version later.
% Set epsilon to a larger value so checking the gradient numerically makes sense
epsilon = 1e-2;
% We pass in the identity matrix as the grouping matrix, putting each
% feature in a group on its own.
groupMatrix = eye(numFeatures);
[cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix);
numgrad = computeNumericalGradient( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix), featureMatrix(:) );
% Uncomment the blow line to display the numerical and analytic gradients side by side
% disp([numgrad grad]);
diff = norm(numgrad-grad)/norm(numgrad+grad);
fprintf('Feature difference (non-topographic): %g\n', diff);
assert(diff < 1e-8, 'Feature difference too large. Check your feature cost function. ');
%% STEP 2c: Implement and test feature cost (topographic)
% Implement sparseCodingFeatureCost in sparseCodingFeatureCost.m and check
% the gradient. This time, we will pass a random grouping matrix in to
% check if your costs and gradients are correct for the topographic
% version.
% Set epsilon to a larger value so checking the gradient numerically makes sense
epsilon = 1e-2;
% This time we pass in a random grouping matrix to check if the grouping is
% correct.
groupMatrix = rand(100, numFeatures);
[cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix);
numgrad = computeNumericalGradient( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix), featureMatrix(:) );
% Uncomment the blow line to display the numerical and analytic gradients side by side
% disp([numgrad grad]);
diff = norm(numgrad-grad)/norm(numgrad+grad);
fprintf('Feature difference (topographic): %g\n', diff);
assert(diff < 1e-8, 'Feature difference too large. Check your feature cost function. ');
%%======================================================================
%% STEP 3: Iterative optimization
% Once you have implemented the cost functions, you can now optimize for
% the objective iteratively. The code to do the iterative optimization
% using mini-batching and good initialization of the features has already
% been included for you.
%
% However, you will still need to derive and fill in the analytic solution
% for optimizing the weight matrix given the features.
% Derive the solution and implement it in the code below, verify the
% gradient as described in the instructions below, and then run the
% iterative optimization.
% Initialize options for minFunc
options.Method = 'lbfgs';
options.display = 'off';
options.verbose = 0;
% Initialize matrices
weightMatrix = rand(visibleSize, numFeatures);
featureMatrix = rand(numFeatures, batchNumPatches);
% Initialize grouping matrix
assert(floor(sqrt(numFeatures)) ^2 == numFeatures, 'numFeatures should be a perfect square');
donutDim = floor(sqrt(numFeatures));
assert(donutDim * donutDim == numFeatures,'donutDim^2 must be equal to numFeatures');
groupMatrix = zeros(numFeatures, donutDim, donutDim);
groupNum = 1; %% 获得拓扑稀疏编码 这段处理不太懂啊！！
for row = 1:donutDim
for col = 1:donutDim
groupMatrix(groupNum, 1:poolDim, 1:poolDim) = 1;
groupNum = groupNum + 1;
groupMatrix = circshift(groupMatrix, [0 0 -1]);
end
groupMatrix = circshift(groupMatrix, [0 -1, 0]);
end
groupMatrix = reshape(groupMatrix, numFeatures, numFeatures);
if isequal(questdlg('Initialize grouping matrix for topographic or non-topographic sparse coding?', 'Topographic/non-topographic?', 'Non-topographic', 'Topographic', 'Non-topographic'), 'Non-topographic')
groupMatrix = eye(numFeatures);
end
% Initial batch
indices = randperm(numPatches);
indices = indices(1:batchNumPatches);
batchPatches = patches(:, indices);
fprintf('%6s%12s%12s%12s%12s\n','Iter', 'fObj','fResidue','fSparsity','fWeight');
for iteration = 1:200 %% 因为要交替优化直到最小化cost function, 所以才这样进行的
error = weightMatrix * featureMatrix - batchPatches;
error = sum(error(:) .^ 2) / batchNumPatches;
fResidue = error;
R = groupMatrix * (featureMatrix .^ 2);
R = sqrt(R + epsilon);
fSparsity = lambda * sum(R(:));
fWeight = gamma * sum(weightMatrix(:) .^ 2);
fprintf(' %4d %10.4f %10.4f %10.4f %10.4f\n', iteration, fResidue+fSparsity+fWeight, fResidue, fSparsity, fWeight) %% 以上这部分可以不用的，只是为了显示最终的
% Select a new batch
indices = randperm(numPatches); %% 重新挑选2000个样本用来进行训练
indices = indices(1:batchNumPatches);
batchPatches = patches(:, indices); %%% 重新挑选的样本
% Reinitialize featureMatrix with respect to the new batch
featureMatrix = weightMatrix' * batchPatches; %% trick 对featureMatrix(s)进行初始化 --技巧 方法
normWM = sum(weightMatrix .^ 2)'; %%%%% 也就是weightMatrix矩阵每列的平方和
featureMatrix = bsxfun(@rdivide, featureMatrix, normWM); %% featureMatrix除以上者
% Optimize for feature matrix
options.maxIter = 20; % 迭代20次，并对featureMatrix进行无约束优化
[featureMatrix, cost] = minFunc( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix), ...
featureMatrix(:), options);
featureMatrix = reshape(featureMatrix, numFeatures, batchNumPatches);
% Optimize for weight matrix
weightMatrix = zeros(visibleSize, numFeatures); %%% 通过直接求导得出对weightMatrix进行优化，这里无需进行梯度迭代或者牛顿法等得出最终的结果
weightMatrix = batchPatches*featureMatrix'/(gamma*batchNumPatches* eye(size(featureMatrix, 1)) + featureMatrix*featureMatrix');
% -------------------- YOUR CODE HERE --------------------
% Instructions:
% Fill in the analytic solution for weightMatrix that minimizes
% the weight cost here.
% Once that is done, use the code provided below to check that your
% closed form solution is correct.
% Once you have verified that your closed form solution is correct,
% you should comment out the checking code before running the
% optimization.
[cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix);
assert(norm(grad) < 1e-12, 'Weight gradient should be close to 0. Check your closed form solution for weightMatrix again.')
error('Weight gradient is okay. Comment out checking code before running optimization.');
% -------------------- YOUR CODE HERE --------------------
% Visualize learned basis
figure(1);
display_network(weightMatrix);
end
sparseCodingWeight.m
function [cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix)
%sparseCodingWeightCost - given the features in featureMatrix,
% computes the cost and gradient with respect to
% the weights, given in weightMatrix
% parameters
% weightMatrix - the weight matrix. weightMatrix(:, c) is the cth basis
% vector.
% featureMatrix - the feature matrix. featureMatrix(:, c) is the features
% for the cth example
% visibleSize - number of pixels in the patches
% numFeatures - number of features
% patches - patches
% gamma - weight decay parameter (on weightMatrix)
% lambda - L1 sparsity weight (on featureMatrix)
% epsilon - L1 sparsity epsilon
% groupMatrix - the grouping matrix. groupMatrix(r, :) indicates the
% features included in the rth group. groupMatrix(r, c)
% is 1 if the cth feature is in the rth group and 0
% otherwise.
if exist('groupMatrix', 'var')
assert(size(groupMatrix, 2) == numFeatures, 'groupMatrix has bad dimension');
else
groupMatrix = eye(numFeatures);
end
numExamples = size(patches, 2);
weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);
featureMatrix = reshape(featureMatrix, numFeatures, numExamples);
% -------------------- YOUR CODE HERE --------------------
% Instructions:
% Write code to compute the cost and gradient with respect to the
% weights given in weightMatrix.
% -------------------- YOUR CODE HERE --------------------
ave_square = sum(sum((weightMatrix * featureMatrix - patches).^2))./numExamples; %计算重构误差
weight_decay = gamma * sum(sum(weightMatrix.^2)); %
cost = ave_square + weight_decay;
grad = (2*weightMatrix*featureMatrix*featureMatrix' - 2 * patches*featureMatrix')./numExamples + 2*gamma*weightMatrix;
grad = grad(:);
endsparseCodingFeatureCost.m
function [cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix)
%sparseCodingFeatureCost - given the weights in weightMatrix,
% computes the cost and gradient with respect to
% the features, given in featureMatrix
% parameters
% weightMatrix - the weight matrix. weightMatrix(:, c) is the cth basis
% vector.
% featureMatrix - the feature matrix. featureMatrix(:, c) is the features
% for the cth example
% visibleSize - number of pixels in the patches
% numFeatures - number of features
% patches - patches
% gamma - weight decay parameter (on weightMatrix)
% lambda - L1 sparsity weight (on featureMatrix)
% epsilon - L1 sparsity epsilon
% groupMatrix - the grouping matrix. groupMatrix(r, :) indicates the
% features included in the rth group. groupMatrix(r, c)
% is 1 if the cth feature is in the rth group and 0
% otherwise.
isTopo = 1;
if exist('groupMatrix', 'var')
assert(size(groupMatrix, 2) == numFeatures, 'groupMatrix has bad dimension');
if(isequal(groupMatrix, eye(numFeatures)))
isTopo = 0;
end
else
groupMatrix = eye(numFeatures);
isTopo = 0;
end
numExamples = size(patches, 2);
weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);
featureMatrix = reshape(featureMatrix, numFeatures, numExamples);
% -------------------- YOUR CODE HERE --------------------
% Instructions:
% Write code to compute the cost and gradient with respect to the
% features given in featureMatrix.
% You may wish to write the non-topographic version, ignoring
% the grouping matrix groupMatrix first, and extend the
% non-topographic version to the topographic version later.
% -------------------- YOUR CODE HERE --------------------
ave_square = sum(sum((weightMatrix * featureMatrix - patches).^2))./numExamples; % 计算重构误差
sparsity = lambda .* sum(sum(sqrt( groupMatrix * (featureMatrix.^2) + epsilon))); %计算系数惩罚项
cost = ave_square + sparsity;
gradResidue = (2* weightMatrix'* weightMatrix*featureMatrix - 2*weightMatrix'*patches)./numExamples; %%+ lambda*featureMatrix./sqrt(featureMatrix.^2+epsilon);
if ~isTopo
gradSparsity = lambda*featureMatrix./sqrt(featureMatrix.^2+epsilon); %%% 不是拓扑的稀疏编码
else
gradSparsity = lambda * groupMatrix'*(groupMatrix *(featureMatrix .^ 2) + epsilon).^(0.5).*featureMatrix; %% 拓扑稀疏编码
end
grad = gradResidue + gradSparsity;
grad = grad(:);
end参考文献：
1：UFLDL教程http://ufldl.stanford.edu/wiki/index.php/UFLDL%E6%95%99%E7%A8%8B
2：http://blog.csdn.net/zouxy09/article/details/24971995/机器学习中的范数规则化之（一）L0、L1与L2范数
3：http://www.cnblogs.com/tornadomeet/archive/2013/04/16/3024292.htmlDeep learning：二十九(Sparse coding练习)
4：http://www.cnblogs.com/tornadomeet/archive/2013/04/13/3018393.htmlDeep learning：二十六(Sparse coding简单理解)
UFLDL教程笔记及练习答案六（稀疏编码与稀疏编码自编码表达）
原文地址：http://blog.csdn.net/lu597203933/article/details/46489647
