UFLDL稀疏编码器练习实现

来源：互联网发布：星际穿越影评知乎编辑：程序博客网时间：2024/05/17 05:55

今天终于做完UFLDL稀疏编码器单元练习，回头看看也没有什么，但调试程序时不知道怎么就那么多问题。归纳起来，有两点值得以后注意：

1，计算KL散度时，作为参考的标量rho要放在上面，而rhohat 向量要放在下面。这本来不是个问题，但使得程序只能正确算出W2，b2梯度值，算不出W1，b1值，通不过梯度验证，浪费好多时间，至今还琢磨不透原因。

2，梯度检验时，逐参数、逐样本的梯度计算是最费时间的。想用全局变量改进一下，但未钻研后面的minFunc函数，怕影响优化，作罢。

3，完成这个作业用时并不长，程序整体运行约15分钟。其中梯度检验约12分钟，优化约2分钟。

下面是我的全部程序。其中minFunc函数未钻研分析，直接套用作业中提供源码。其它函数基本都依作业内容翻译或者改写，并完成了要求的功能代码。

%% CS294A/CS294W Programming Assignment Starter Code

% 说明
% ------------
%
% 这个文件包含帮助你开始编程作业的代码。你得完成sampleIMAGES.m
% sparseAutoencoderCost.m 以及 computeNumericalGradient.m中代码
% 就完成作业而言，不用改变文件代码
%
%%======================================================================
%% STEP 0: 这里提供相关参数值，使你的稀疏自编码器能有好的滤波器
% 不用改变下面参数

visibleSize       = 8*8;          % 输入单元数量：相应于8x8图像块
hiddenSize      = 25;            % 隐藏单元数量：25
sparsityParam = 0.01;        % 想要的隐藏单元平均活动度
lambda            = 0.0001;     % 权值衰减参数
beta                 = 3;              % 稀疏权值惩罚项

%%======================================================================
%% STEP 1: 实现 sampleIMAGES
%
% 实现 sampleIMAGES之后，display_network命令应该显示
% 来自数据集容量为200的随机样本块

patches = sampleImagesWEI;
img200 = patches(:,randi(size(patches,2),200,1));
display_networkWEI(img200, [], [], 20);

% 运行程序得到稀疏自编码器（W1，W2，b1，b2）一组随机初始化值
theta = initializeParametersWEI(hiddenSize, visibleSize);

%%======================================================================
%% STEP 2: 实现 sparseAutoencoderCost
%
% 可以一次实现代价函数中所有成分（平方误差代价项，权值衰减项，稀疏性惩罚项），
% 但一步步做，每步后运行梯度检验（见STEP 3）可能会容易些。建议按下述步骤实现
% sparseAutoencoderCost:
%
% (a) 在神经网络中实现前向传导，并实现代价函数的平方误差项
% 实现后向传导计算导数，然后取参数lambda=beta=0
% 运行梯度检验，证实相应于平方误差代价项的计算是正确的
%
% (b) 加入权值衰减项（代价函数和导数计算中都加入），再用梯度检验证实正确性
%
% (c) 加入稀疏性惩罚项，再用梯度检验证实正确性
%
% 调试代码时随便修改训练设置
% (例如，缩减训练集或隐藏单元数能加快代码运行，beta 与 lambda置0有助于调试)
% 不过，在最终提交的可视化权值中，请使用上面STEP 0中提供参数

[~, grad] = sparseAutoencoderCostWEI(theta, visibleSize, hiddenSize, ...
lambda, sparsityParam, beta, patches);

%%======================================================================
%% STEP 3: 梯度检验
%
% 提示：调试代码时，在小模型和小训练集上执行梯度检验能加快速度
% （例如，仅用10个训练样本和1-2个隐藏单元）

% 首先，确信数值方法计算的梯度在简单函数上是正确的
% 在已经实现 computeNumericalGradient.m 之后，执行如下命令:
checkNumericalGradientWEI();

% 现在就可以用它检验稀疏自编码器的代价函数和导数计算
tic
numgrad = computeNumericalGradientWEI( @(x)   ...
             sparseAutoencoderCostWEI(x, visibleSize, hiddenSize, ...
                lambda, sparsityParam, beta, patches), theta);

% 用下面函数可视化比较两个梯度：方便起见，画出两个梯度只差的曲线
% disp([numgrad grad]);
plot(numgrad-grad)

% 比较数值计算方法获得梯度与后向传导方法获得梯度：完全一致，误差<1e-10
diff = norm(numgrad-grad)/norm(numgrad+grad);
disp(diff); % 应该很小。我们的实现中，这些值通常小于 1e-9.
% 如果做到这点，祝贺 !!!
toc

%%======================================================================
%% STEP 4: 在确信你的sparseAutoencoderCost实现正确后，
% 即可开始用 minFunc(L-BFGS) 训练稀疏编码器
tic
% 随机初始化参数
theta = initializeParameters(hiddenSize, visibleSize);

% 用 minFunc 最小化函数
addpath minFunc/
options.Method = 'lbfgs'; % 此处用 L-BFGS 优化代价函数。minFunc通常需要
                          % 有两个输出的函数指针：函数值与梯度
                          % sparseAutoencoderCost.m 满足条件
options.maxIter = 400;   % 运行 L-BFGS 最大循环数值
options.display = 'on';

[opttheta, cost] = minFunc( @(p) sparseAutoencoderCostWEI(p, ...
                                   visibleSize, hiddenSize, ...
                                   lambda, sparsityParam, ...
                                   beta, patches), ...
                              theta, options);
toc
%%======================================================================
%% STEP 5: 可视化

W1 = reshape(opttheta(1:hiddenSize*visibleSize), hiddenSize, visibleSize);
display_networkWEI(W1', 12);

print -djpeg weights.jpg % 图像保存到文件中

function patches = sampleImagesWEI()
% sampleIMAGES，返回 10000 个训练用图像样本块，64x10000 矩阵形式
% 图像样本块是从所提供图像中随机剪切出的 8x8 块

load IMAGES;                % 装入所提供图像
patchsize = 8;             % 定义图像快大小（8x8）
numpatches = 10000;         % 图像块数量（10000）

% 所有图像块初始化为0，下面编写代码从图像中剪切图像块填充
patches = zeros(patchsize*patchsize, numpatches);

%% ---------- YOUR CODE HERE --------------------------------------
% 说明：用来自 IMAGES 的数据填充称为 patches 的变量
%
% IMAGES 是3D数组，包含10幅图像
% 举例：IMAGES(:,:,6) 是第6幅图像，大小为 512x512
% 可用 imagesc(IMAGES(:,:,6)), colormap gray;命令看看这幅图像
% (这些图像的对比度有些暗，因为被用“白化”方法预处理过，细节见讲义）
% 再举一例：IMAGES(21:30,21:30,1) 是一个 10x10 图像块，对应第一幅图像中
% （21，21）到（30，30）区域图像

% randi(IMax,1,N)生成 1xN 随机向量，元素在 [1 IMAX]中均匀采样，有重复

image_size = size(IMAGES);
rowStarter = randi(image_size(1) - patchsize + 1, 1, numpatches);
colStarter = randi(image_size(2) - patchsize + 1, 1, numpatches);
imgSelected = randi(image_size(3), 1, numpatches);

for k = 1 : numpatches
    patches(:, k) = reshape(          ...
        IMAGES( rowStarter(k) : rowStarter(k) + patchsize - 1, ...
                colStarter(k) : colStarter(k) + patchsize - 1, ...
                imgSelected(k) ),      ...
        1, patchsize * patchsize);
end

%% ---------------------------------------------------------------
% 为让自编码器好好工作，需要归一化数据
% 这里，因为网络输出界定在 [0,1] 之间（sigmoid激活函数决定）
% 故得确保像素值也在 [0,1] 之间
patches = normalizeData(patches);

end

%% ---------------------------------------------------------------
function patches = normalizeData(patches)

% 因为在输出层使用sigmoid激活函数，故将数据压入 [0.1, 0.9]

% 去除直流成分（图像均值）
patches = bsxfun(@minus, patches, mean(patches));

% 截尾为3倍的标准离差，调整为 -1 到 1
pstd = 3 * std(patches(:));
patches = max(min(patches, pstd), -pstd) / pstd;

% 从 [-1,1] 再调整为 [0.1,0.9]
patches = (patches + 1) * 0.4 + 0.1;

end

function theta = initializeParametersWEI(hiddenSize, visibleSize)

%% 基于层的宽度随机初始化参数

% 权值W1、W2初始值均匀分布在[-r, r]，偏执值 b1、b2初始值取0
r = sqrt(6) / sqrt(hiddenSize+visibleSize+1);
W1 = rand(hiddenSize, visibleSize) * 2 * r - r;
W2 = rand(visibleSize, hiddenSize) * 2 * r - r;

b1 = zeros(hiddenSize, 1);
b2 = zeros(visibleSize, 1);

% 权值与偏置值转换为向量形式，所有参数按minFunc要求，展开合并为一个向量
theta = [W1(:) ; W2(:) ; b1(:) ; b2(:)];

end

function [cost,grad]=sparseAutoencoderCostWEI(theta, visSize, ...
hidSize, lambda, sparsityParam, beta, data)

% 计算稀疏自编码器的代价函数值与梯度值
sampleNum=size(data,2);

% 从theta向量恢复W1，W2矩阵，b1，b2向量
W1 = reshape(theta(1:hidSize*visSize), hidSize, visSize);
W2 = reshape(theta(hidSize*visSize+1:2*hidSize*visSize), ...
visSize, hidSize);
b1 = theta(2*hidSize*visSize+1:2*hidSize*visSize+hidSize);
b2 = theta(2*hidSize*visSize+hidSize+1:end);

% 计算代价函数：隐层与输出层活动度
z2 = W1*data + repmat(b1,1,sampleNum);
a2 = sigmoid(z2);
z3 = W2*a2 + repmat(b2,1,sampleNum);
a3 = sigmoid(z3);
rho = sum(a2,2)./sampleNum;
cost = 1/2/sampleNum * sum(sum((a3-data).^2)) + ...
lambda/2 * (sum(sum(W1.^2))+sum(sum(W2.^2))) + ...
beta * sum(KLDivergence(sparsityParam,rho));

% 计算输出层与隐层残差，根据残差计算梯度值
sparseX = beta*(-sparsityParam./rho+(1-sparsityParam)./(1-rho));
delta3 = -(data-a3).*a3.*(1-a3);
delta2 = (W2'*delta3+repmat(sparseX,1,sampleNum)).*a2.*(1-a2);

W1grad = delta2*data'/sampleNum + lambda*W1 ;
W2grad = delta3*a2'/sampleNum + lambda*W2;
b1grad = sum(delta2,2)/sampleNum;
b2grad = sum(delta3,2)/sampleNum;

grad = [W1grad(:); W2grad(:); b1grad(:); b2grad(:)];
end

function sigm=sigmoid(x)
% sigmoid 传导函数
sigm = 1 ./ (1 + exp(-x));
end

function entropyX=KLDivergence(rho,rhoHat)

% 计算KL(rho||rhohat)
entropyX = rho.*log(rho./rhoHat)+(1-rho).*log((1-rho)./(1-rhoHat));
end

function numgrad = computeNumericalGradientWEI(J, theta)
% numgrad = computeNumericalGradient(J, theta)
% theta: 参数向量
% J: 函数，输出实数值
% y = J(theta) 返回函数在theta处的值

% Initialize numgrad with zeros
numgrad = zeros(size(theta));

%% ---------- YOUR CODE HERE -----------------------------------
% 说明
% 实现数值梯度检验，结果返回在numgrad (参考讲义 2.3节）
% 编写代码使 numgrad(i) 是J相对于第i个输入项在theta处的偏导数数值方法近似值
% 即 numgrad(i) 应为J相对于theta(i)的偏导数近似值
%
% 提示：一次计算一个梯度

% 按公式计算
EPSILON = 0.0001;
prmtNum = length(theta);
sqrDiag = eye(prmtNum);
for i = 1 : prmtNum
delta = sqrDiag(:,i)*EPSILON;
numgrad(i) = (J(theta+delta)-J(theta-delta))/(2*EPSILON);
end

%% -------------------------------------------------------------
end

function [] = checkNumericalGradientWEI()

% 对比函数h(x1,x2)=x1.^2+3*x1*x2数值方法梯度值与理论梯度值

% 选函数值与梯度比较点：x = [4; 10]
x = [4; 10];
[~,grad] = simpleQuadraticFunction(x);

% 编写代码用数值方法计算 simpleQuadraticFunction 在 x 处的值与梯度
numgrad = computeNumericalGradientWEI(@simpleQuadraticFunction, x);

% 可视化两个梯度，应该很像
disp([numgrad grad]);
fprintf('(Left - Numerical Gradient, Right - Analytical Gradient)\n\n');

% 评估两个解差异的范数。如果 EPSILON = 0.0001，则 diff < 2.1452e-12
diff = norm(numgrad-grad)/norm(numgrad+grad);
disp(diff);
fprintf('Norm of the difference should be less than 1e-9 \n\n');

end

function [value,grad] = simpleQuadraticFunction(x)
% 输入2D向量x=（x1，x2）
% 输出
% value: h(x1, x2) = x1^2 + 3*x1*x2
% grad: 2x1 向量，h相对于 x1 与 x2偏导数
% 注意，用 @simpleQuadraticFunction(x)方法调用时，假定只用到第一个返回值

value = x(1)^2 + 3*x(1)*x(2);
grad = zeros(2, 1);
grad(1) = 2*x(1) + 3*x(2);
grad(2) = 3*x(1);

end

function [h, array] = display_networkWEI(A, opt_normalize, ...
opt_graycolor, cols, opt_colmajor)
% 这个函数可视化矩阵A的滤波器。A的每一列是一个滤波器，变形为一个方阵
% 可视化在图像面板每一个格子中
% 所有其它参数都是可选的，通常不用管它们
% opt_normalize: 要不要归一化滤波器，以便全部有统一的对比度，默认为真
% opt_graycolor: 是不是用灰度作为热图，默认是
% cols: 显示格子有多少列？默认值为A列数的平方根
% opt_colmajor: 以列为主或以行为主习惯之间的交互，默认为假

% 原来程序在显示200幅图像时，显示为12x17（=204）排列，最后4个为空，不理想

warning off all

if ~exist('opt_normalize', 'var') || isempty(opt_normalize)
opt_normalize= true;
end

if ~exist('opt_graycolor', 'var') || isempty(opt_graycolor)
opt_graycolor= true;
end

if ~exist('opt_colmajor', 'var') || isempty(opt_colmajor)
opt_colmajor = false;
end

% 重新调节元素值
A = A - mean(A(:));

if opt_graycolor, colormap(gray); end

% 计算行数m，列数n
[L M]=size(A);
sz=sqrt(L);
buf=1;
if ~exist('cols', 'var')
    if floor(sqrt(M))^2 ~= M
        n=ceil(sqrt(M));
        while mod(M, n)~=0 && n<1.2*sqrt(M), n=n+1; end
        m=ceil(M/n);
    else
        n=sqrt(M);
        m=n;
    end
else
    n = cols;
    m = ceil(M/n);
end

array=-ones(buf+m*(sz+buf),buf+n*(sz+buf));

if ~opt_graycolor
array = 0.1.* array;
end

if ~opt_colmajor
    k=1;
    for i=1:m
        for j=1:n
            if k>M,
                continue;
            end
            clim=max(abs(A(:,k)));
            if opt_normalize
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/clim;
            else
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/max(abs(A(:)));
            end
            k=k+1;
        end
    end
else
    k=1;
    for j=1:n
        for i=1:m
            if k>M,
                continue;
            end
            clim=max(abs(A(:,k)));
            if opt_normalize
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/clim;
            else
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz);
            end
            k=k+1;
        end
    end
end

if opt_graycolor
h=imagesc(array,'EraseMode','none',[-1 1]);
else
h=imagesc(array,'EraseMode','none',[-1 1]);
end
axis image off

drawnow;

warning on all

0 0