UFLDL作业记录

UFLDL-Sparse Autoencoder

稀疏自编码的作业，在这里。

1. 实现sampleIMAGES

IMAGES.mat是一个512x512x10的10张已将白化过的图片。

sampleIMAGES，故名思议，就是采样IMAGES。

用

imagesc(IMAGES(:,:,4)), colormap gray

可以看到图像的灰度的样子,比如第四张图片：



我们想要采的样本是10000个8x8的样本，因为sampleIMAGES里已经这么初始化了。

patchsize = 8;  % we'll use 8x8 patches
numpatches = 10000;
patches = zeros(patchsize*patchsize, numpatches);

patches被初始化程成了一个64x10000的矩阵。

Idea很简单：

随机从10张图片里选一张，随机从x坐标里选8个像素，随机从y坐标里选8个像素。

其实我们只要定下10000个随机的点，作为8x8的样本的左上角就可以了。

设，左上角的点的坐标是(a,b,c)，显然(a+8-1,b+8-1,c)要不大于(512,512,10)

先选x坐标的值：

a = randi(512-8+1,1,10000);

在选y坐标的值：

b = randi(512-8+1,1,10000);

最后选图片到底是第几张：

c = randi(10,1,10000);

于是，10000张8x8的图像，其实就是：

d = IMAGES(a:a+8-1,b:b+8-1,c)

所以patches就是：

for i = 1:10000
pathces(:,i) = reshape(d(:,:,i),1,64);
end

综合起来，就是：

tic
image_size = size(IMAGES);
a = randi(image_size(1) - patchsize + 1, 1, numpatches);
b = randi(image_size(2) - patchsize + 1, 1, numpatches);
c = randi(image_size(3), 1, numpatches);
d = IMAGES(a:a+8-1, b:b+8-1, c)
for num = 1 : numpatches
patches(:,num)=reshape(d(:,:,num),1,patchsize*patchsize);
end
toc

2. 构建Sparse autoencoder objective

可以看到，train.m里，已经设定了hiddenSize = 25，所以这个网络结构其实是：64x25x64，隐层有25个单元。

先按正常的神经网络的训练顺序来

2.1 前向传播

Θ已经在initializeParameters里定义好了，所以 W~1~(25x64),W~2~(64x25)和b~1~(25x1),b~2~(64x1) 已经有了，可以直接开始算前向传播：

B_1 = repmat(b_1,1,10000);%截距
Z_2 = (W_1*patches) + B_1;%输入
A_2 = sigmoid(Z_2);%隐层输出
B_2 = repmat(b_2,1,10000);%截距
Z_3 = (W_2*A_2) + B_2;%输入
A_3 = sigmoid(Z_3);%结果

2.2 稀疏性保证

接着要计算ρ^,号称是隐层的平均输出，文中给的公式是：

ρ^j=1m∑mi=1[a(2)j(x(i))]

它说，这里取的平均，是在测试集上的平均(averaged over the training set),我觉得很奇怪。

我们知道A_2是一个25x10000的矩阵，那这个m，到底是25还是10000？我怎么感觉都是25，不是说是隐层的平均输出么？

我们想让ρ^尽量的小，最好是某个挺小的值ρ,这个ρ也叫稀疏性参数

为了完成上述要求，我们要把ρ^和ρ的相对熵算一下：

KL(ρ||ρ^j)=ρlogρρ^j+(1−ρ)log1−ρ1−ρ^j

把下面这个东西放进代价函数里，来实现稀疏性的控制：

∑s2j=1KL(ρ||ρ^j)

文中说，说这里的S~2~是隐层节点数，也就是25。这里的ρ又是个标量，显然上面的那个ρ^里的m是10000了。

那么，代码这么写：

mean_A_2 = sum(A_2,2)./10000;
rho = sparsityParam
KL = sum(rho.*log(rho./mean_A_2)+(1-rho).*log((1-v)/(1.-mean_A_2));

2.3 代价函数

代价函数J：

Jsparse(W,b)=J(W,b)+β∑s2j=1KL(ρ||ρ^j)

其中，J(W,b),按照反向传播来算是：

J(W,b)=[1m∑mi=1(12∥∥hW,b(x(i))−y(i)∥∥2)]+λ2∑nl−1l=1∑sli=1∑sl+1j=1(W(l)ji)2

所以，

J=(1/2*10000)*sum(sum((A_3-patches).^2))+lambda/2*(sum(sum(W_1.^2))+sum(sum(W_2.^2)))+beta*KL;

2.4梯度下降

daoshu_Z_3 = sigmoid(A_3).*(1-sigmoid(A_3);
delta_3 = -(patches-A_3).*daoshu_Z_3;

对于delta_2,文中说：

δ(2)i=((∑s2j=1W(2)jiδ(3)j)+β(−ρρ^i+1−ρ1−ρ^i))f′(z(2)i)

这么写比较省事。所以：

sparity_penalty = (-rho./mean_A_2+(1-rho)./(1-mean_A_2));
daoshu_Z_2 = sigmoid(A_2).*(1-sigmoid(A_2);
delta_2=(W_2'*delta_3+repmat(beta*sparity_penalty,1,mm)).*daoshu_Z_2;

则，W偏导，b偏导为：

∇W(l)J(W,b;x,y)=δ(l+1)(a(l))T

∇b(l)J(W,b;x,y)=δ(l+1)

所以W,和b的偏导代码是：

W1_pd = delta_2 * pathces';
W2_pd = delta_3 * A_2';
b1_pd = sum(delta_2,2);
b2_pd = sum(delta_3,2);

W，b的梯度：

ΔW(l):=ΔW(l)+∇W(l)J(W,b;x,y)

Δb(l):=Δb(l)+∇b(l)J(W,b;x,y)

W1grad = W1grad + W1_pd;
W2grad = W2grad + W2_pd;
b1grad = b1grad + b1_pd;
b2grad = b2grad + b2_pd;

最终更新权重：

b1grad=b1grad/10000;
b2grad=b2grad/10000;
W1grad=W1grad/10000+lambda*W_1;
W2grad=W2grad/10000+lambda*W_2;

至此，代价函数就写完了。

3.检验梯度

最后，还要保证梯度算的没错。。。。

文中提到了这样的检测方法：

g(θ)≈J(θ+EPSILON)−J(θ−EPSILON)2×EPSILON.

那么，我们就要对代价函数做上述的检测：

eps = 1e-4;
m = size(theta,1);
for i=1:m
theta_p = theta;
theta_p(i,:) = theta_p(i,:) + eps;
theta_m = theta;
theta_m(i,:) = theta_m(i,:) - eps;
numgrad(i,:) = (J(theta_p) - J(theta_m))/(eps*2.0);
end

其实有向量化的方法：

eps = 1e-4;
m = size(theta,1);
E = eye(m);
for i = 1:m
delta = E(:,i)*eps;
numgrad(i) = (J(theta+delta)-J(theta-delta))/(eps*2.0);
end

m一大，eye就超出了matlab的限制了，所以，我用的上一个。

这一步句，巨慢无比，因为要进行mx2次的J函数的运算。。。。

4.结果