Unsupervised Feature Learning and Deep Learning(UFLDL) Exercise 代码实现

最近开始了解深度学习的东西，看到国内研究者翻译的斯坦福大学的深度学习的历程，觉得很不错，就按照学习例程学习，一步一步完成后面的练习。这个练习安排的很合理，得按顺序完成，跳跃着做可能会遇到一定的麻烦。我在完成中也遇到不少的麻烦，不过还是慢慢克服了，越往后程序对计算机的配置要求越高，需要的内存也越大。下面是我完成这些例程中一些关键性的代码，觉得有用的可以参考一下，建议大家先自己动手写，如果实在是遇到困难了，可以参考一下我写的代码，由于本人编程能力也很一般，代码优化做得也很差，只供参考而已。

Exercise1: Sparase AutoEncoder

这个练习比较关键，基本上后面许多练习都依赖于这个练习的代码，所以首先得完成好这个练习。这个练习需要补充3个地方的代码。我这里只列出其中最关键的一个地方的代码，也就是sparseAutoencoderCost.m这个文件。
函数的目的是首先代价的计算以及梯度的计算，具体的实现如下：

A1 = data;
Y = data;
z2 = W1 * data + repmat(b1,1,m);
A2 = sigmoid(z2);
z3 = W2 * A2 + repmat(b2,1,m);
A3 = sigmoid(z3);

rho_hat = mean(A2,2); % 隐层每个神经元的平均激活度 hiddenSize*1

%backpropagation
for i = 1:m

a1 = A1(:,i);

a2 = A2(:,i);

a3 = A3(:,i);

y = Y(:,i);                 % 原始输入向量

cost = cost + 1/2 * norm(y - a3, 2)^2;   % 计算代价

delta_3 = -(y - a3) .* (a3 .* (1 - a3));

delta_2 = (W2'* delta_3 + beta *(-sparsityParam ./rho_hat + (1 - sparsityParam) ./(1 - rho_hat))) .* (a2 .* (1 - a2));

W1grad = W1grad + delta_2 * a1';
b1grad = b1grad + delta_2;
W2grad = W2grad + delta_3 * a2';
b2grad = b2grad + delta_3;
end

% 计算梯度
W1grad = 1/m * W1grad + lambda * W1;
W2grad = 1/m * W2grad + lambda * W2;
b1grad = 1/m * b1grad;
b2grad = 1/m * b2grad;

%计算KL
KL = 0;
for i = 1:hiddenSize
KL = KL + sparsityParam * log(sparsityParam/rho_hat(i)) + (1 - sparsityParam) * log((1 - sparsityParam)/(1 - rho_hat(i)));
end

%计算代价
cost = 1/m * cost + lambda/2 * (sum(sum(W1.^2)) + sum(sum(W2.^2))) + beta * KL;</span>

注意：自己实现的时候一定要按照练习提示一步一步来，基本上可以搞定的，还有，我这个实现其实不是最开始的代码，这里面参杂了一部分练习2的代码，不过练习2基本上是完全优化了上面的代码，去掉了所有影响运行效率的for循环，但是带来的问题是程序的内存开销也相应增大了很多，所以这点得注意，不然会报out of memory的错误。

Exercise2:Vectorization

这个练习我个人觉得很重要，它交给我们如何将for循环尽量改写为向量或矩阵运算的方式，对程序运行速度有大幅的提升。这里的代码是对Exercise1的重新实现：

z2 = W1 * data + repmat(b1,1,m);
a2 = sigmoid(z2);
z3 = W2 * a2 + repmat(b2,1,m);
a3 = sigmoid(z3);

rho_hat = mean(a2,2); % 隐层每个神经元的平均激活度 hiddenSize*1
a1 = data;
y = data;

cost = 1/2*sum(sum((y - a3).^2));
delta_3 = -(y - a3) .* (a3 .* (1 - a3));
sparsity_data = (-sparsityParam ./rho_hat + (1 - sparsityParam) ./(1 - rho_hat));
delta_2 = (W2'* delta_3 + beta * repmat(sparsity_data,1,m)) .* (a2 .* (1 - a2));
W1grad = delta_2 * a1';
b1grad = sum(delta_2,2);
W2grad = delta_3 * a2';
b2grad = sum(delta_3,2);

W1grad = 1/m * W1grad + lambda * W1;
W2grad = 1/m * W2grad + lambda * W2;
b1grad = 1/m * b1grad;
b2grad = 1/m * b2grad;

% 计算KL
KL = sparsityParam * log(sparsityParam./rho_hat) + (1 - sparsityParam) * log((1 - sparsityParam)./(1 - rho_hat));
KL = sum(KL);

cost = 1/m * cost + lambda/2 * (sum(sum(W1.^2)) + sum(sum(W2.^2))) + beta * KL;

可以看到，没有使用for循环，同样实现了上述的功能。

Exercise3:PCA and Whitening

这个例子和上面的关系不大，其实也不是很难，对PCA有所了解的话实现上该练习完全不在话下。这里面有点需要注意的是对2D图像进行PCA压缩重构的时候，需要用到U和V，而不是只用U就可以，大家可以简单的检查一下矩阵的维度就很容易发现这点。所以下面就只列出重构图像的代码：

k = 0; % Set k accordingly
lambda = 0;
for i = 1:size(x,1)
lambda = lambda + S(i,i);
if lambda/trace(S) >= 0.96
k = i;
break;
end
end

xHat = U(:,1:k) * V(:,1:k)' * x;

Exercise4:Softmax Regression

softmax 回归是Logistic 回归的推广，也是机器学习中常用的分类器。这个实现起来不是很难，但是我们刚刚完成了Exercise2,就得用向量话的方式实现，活学活用嘛，这样也可以逐步锻炼自己的编程习惯。

x = data;
k = numClasses;
m = numCases;
out = exp(theta * x);
S = sum(exp(theta * x),1);

thetagrad = (groundTruth - out ./ repmat(S,k,1)) * x';
thetagrad = -1/m * thetagrad + lambda * theta;

cost = sum(sum(groundTruth .* log(out./repmat(S,k,1))));
cost = -1/m * cost + lambda/2 * sum(sum(theta.^2));

大家看到，是不是和简洁，而且速度很快。所以要多培养哦！

Exercise5:Self-Taught Learning

这个练习基本上是对前面练习结果的一个重用，没什么难度，所以说不能跳着做，得按顺序一步步完成比较好。

Exercise6:Implement Deep Learning Networks for Digital Classification

这个练习可以算的上是自己的第一个真正意义上的深度学习程序，采用了两个自动编码器隐层和一个softmax预测层，其中难点在于实现有监督的调优过程。下面是stackedAECost.m的实现。

k = numClasses;

W1 = stack{1}.w;
b1 = stack{1}.b;
W2 = stack{2}.w;
b2 = stack{2}.b;

output_1 = sigmoid(W1 * data + repmat(b1,1,M));     % 计算第一隐层特征输出
output_2 = sigmoid(W2 * output_1 + repmat(b2,1,M)); % 计算第二隐层特征输出

out = exp(softmaxTheta * output_2);                        % 计算分类器输出
S = sum(exp(softmaxTheta * output_2),1);

% 计算softmax的梯度
softmaxThetaGrad = (groundTruth - out ./ repmat(S,k,1)) * output_2';
softmaxThetaGrad = -1/M * softmaxThetaGrad + lambda * softmaxTheta;

cost = sum(sum(groundTruth .* log(out./repmat(S,k,1))));
cost = -1/M * cost + lambda/2 * sum(sum(softmaxTheta.^2));

% 计算autoEncoder的梯度
delta_4 = -(groundTruth - out ./ repmat(S,k,1));
delta_3 = (softmaxTheta' * delta_4).* (output_2 .* (1 - output_2));
delta_2 = (W2'* delta_3) .* (output_1 .* (1 - output_1));
% delta_1 = (W1' * delta_2) .* (data .* (1 - data));

stackgrad{1}.w = 1/M * delta_2 * data';
stackgrad{1}.b = 1/M * sum(delta_2,2);
stackgrad{2}.w = 1/M * delta_3 * output_1';
stackgrad{2}.b = 1/M * sum(delta_3,2);

注意，由于这个程序使用了两个自动编码器层，所以训练时间是比较长的（配置高的电脑大概也得30min，运行内存最好有8G，不然也可能会报out of memory 的错误。

Exercise7:Linear Decoders with Autoencoders

这个练习和上一个练习基本上是相同的，就是把上个练习中第二个自动编码器的激励函数从sigmoid 换成了线性的，其他的没什么大的变化。

Exercise8:Convolution and Pooling

卷积神经网络在处理自然图像的时候有很大的优势，图像可以直接输入网络而不需要做向量化拉直。这个练习需要用到上个练习中的结果，练习中有提示，所以在做这个练习的时候最好先完成上面的练习。这里面需要实现的是cnnConvolve.m这个文件。注意我在其中有些没有按照他的提示做，但是效果是相同的。

% 计算特征卷积模板
AllFeatures = (W * ZCAWhite)';   % 400 * 192;
convFeatures = reshape(AllFeatures,8,8,3,numFeatures);

% --------------------------------------------------------

convolvedFeatures = zeros(numFeatures, numImages, imageDim - patchDim + 1, imageDim - patchDim + 1);
for imageNum = 1:numImages
for featureNum = 1:numFeatures

% convolution of image with feature matrix for each channel
convolvedImage = zeros(imageDim - patchDim + 1, imageDim - patchDim + 1);
for channel = 1:imageChannels

% Obtain the feature (patchDim x patchDim) needed during the convolution
% ---- YOUR CODE HERE ----
%       feature = zeros(8,8); % You should replace this

feature = squeeze(convFeatures(:,:,channel,featureNum));

% ------------------------

% Flip the feature matrix because of the definition of convolution, as explained later
%feature = flipud(fliplr(squeeze(feature)));
feature = rot90(feature,2);

% Obtain the image
im = squeeze(images(:, :, channel, imageNum));

% Convolve "feature" with "im", adding the result to convolvedImage
% be sure to do a 'valid' convolution
% ---- YOUR CODE HERE ----
convolvedImage = convolvedImage + conv2(im,feature,'valid');

% ------------------------

end

% Subtract the bias unit (correcting for the mean subtraction as well)
% Then, apply the sigmoid function to get the hidden activation
% ---- YOUR CODE HERE ----
bias = b - W*ZCAWhite*meanPatch;
convolvedImage = sigmoid(convolvedImage + repmat(bias(featureNum),imageDim - patchDim + 1, imageDim - patchDim + 1) );

% ------------------------

% The convolved feature is the sum of the convolved values for all channels
convolvedFeatures(featureNum, imageNum, :, :) = convolvedImage;
end
end

end

function y = sigmoid(x)
y = 1 ./ (1 + exp(-x));
end

以上给出的代码是我个人认为比较关键的地方，当然要想让程序跑起来，还得实现其他的地方。所以大家如果在完成这些练习遇到困难实在不知道如何解决的时候可以参考我上面给出的关键代码。