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 前向传播

$\Theta$已经在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 稀疏性保证

接着要计算$\hat\rho$,号称是隐层的平均输出，文中给的公式是：
$\hat\rho_j = \frac{1}{m} \sum_{i=1}^m \left[ a^{(2)}_j(x^{(i)}) \right]$
它说，这里取的平均，是在测试集上的平均(averaged over the training set),我觉得很奇怪。

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

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

为了完成上述要求，我们要把$\hat\rho$和$\rho$的相对熵算一下：
${\rm KL}(\rho || \hat\rho_j) = \rho \log \frac{\rho}{\hat\rho_j} + (1-\rho) \log \frac{1-\rho}{1-\hat\rho_j}$
把下面这个东西放进代价函数里，来实现稀疏性的控制：
$\sum_{j=1}^{s_2} {\rm KL}(\rho || \hat\rho_j)$
文中说，说这里的S~2~是隐层节点数，也就是25。这里的$\rho$又是个标量，显然上面的那个$\hat\rho$里的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：
$J_{\rm sparse}(W,b) = J(W,b) + \beta \sum_{j=1}^{s_2} {\rm KL}(\rho || \hat\rho_j)$
其中，$J(W,b)$,按照反向传播来算是：

$J(W,b)= \left[ \frac{1}{m} \sum_{i=1}^m \left( \frac{1}{2} \left\| h_{W,b}(x^{(i)}) - y^{(i)} \right\|^2 \right) \right]+ \frac{\lambda}{2} \sum_{l=1}^{n_l-1} \; \sum_{i=1}^{s_l} \; \sum_{j=1}^{s_{l+1}} \left( W^{(l)}_{ji} \right)^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,文中说：
$\delta^{(2)}_i = \left( \left( \sum_{j=1}^{s_{2}} W^{(2)}_{ji} \delta^{(3)}_j \right) + \beta \left( - \frac{\rho}{\hat\rho_i} + \frac{1-\rho}{1-\hat\rho_i} \right) \right) 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偏导为：
$\nabla_{W^{(l)}} J(W,b;x,y) = \delta^{(l+1)} (a^{(l)})^T$

$\nabla_{b^{(l)}} J(W,b;x,y) = \delta^{(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的梯度：

$\Delta W^{(l)} := \Delta W^{(l)} + \nabla_{W^{(l)}} J(W,b;x,y)$

$\Delta b^{(l)} := \Delta b^{(l)} + \nabla_{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.检验梯度

最后，还要保证梯度算的没错。。。。
文中提到了这样的检测方法：
$\begin{align} g(\theta) \approx \frac{J(\theta+{\rm EPSILON}) - J(\theta-{\rm EPSILON})}{2 \times {\rm EPSILON}}. \end{align}$
那么，我们就要对代价函数做上述的检测：

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.结果