UFLDL-Sparse Autoencoder
稀疏自编码的作业，在这里。
1. 实现sampleIMAGES
IMAGES.mat是一个512x512x10的10张已将白化过的图片。
sampleIMAGES，故名思议，就是采样IMAGES。
用
imagesc(IMAGES(:,:,)), colormap gray
可以看到图像的灰度的样子,比如第四张图片：
我们想要采的样本是10000个8x8的样本，因为sampleIMAGES里已经这么初始化了。
patchsize = ; % we'll use 8x8 patches
numpatches = ;
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(-+,,);
在选y坐标的值：
b = randi(-+,,);
最后选图片到底是第几张：
c = randi(,,);
于是，10000张8x8的图像，其实就是：
d = IMAGES(a:a+-,b:b+-,c)
所以patches就是：
for i = :
pathces(:,i) = reshape(d(:,:,i),,);
end
综合起来，就是：
tic
image_size = size(IMAGES);
a = randi(image_size() - patchsize + , , numpatches);
b = randi(image_size() - patchsize + , , numpatches);
c = randi(image_size(), , numpatches);
d = IMAGES(a:a+-, b:b+-, c)
for num = : numpatches
patches(:,num)=reshape(d(:,:,num),,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,,);%截距
Z_2 = (W_1*patches) + B_1;%输入
A_2 = sigmoid(Z_2);%隐层输出
B_2 = repmat(b_2,,);%截距
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,)./;
rho = sparsityParam
KL = sum(rho.*log(rho./mean_A_2)+(-rho).*log((-v)/(-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=(/*)*sum(sum((A_3-patches).^))+lambda/*(sum(sum(W_1.^))+sum(sum(W_2.^)))+beta*KL;
2.4梯度下降
daoshu_Z_3 = sigmoid(A_3).*(-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+(-rho)./(-mean_A_2));
daoshu_Z_2 = sigmoid(A_2).*(-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,);
b2_pd = sum(delta_3,);
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/;
b2grad=b2grad/;
W1grad=W1grad/+lambda*W_1;
W2grad=W2grad/+lambda*W_2;
至此，代价函数就写完了。
3.检验梯度
最后，还要保证梯度算的没错。。。。
文中提到了这样的检测方法：
g(θ)≈J(θ+EPSILON)−J(θ−EPSILON)2×EPSILON.
那么，我们就要对代价函数做上述的检测：
eps = ;
m = size(theta,);
for i=: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*);
end
其实有向量化的方法：
eps = ;
m = size(theta,);
E = eye(m);
for i = :m
delta = E(:,i)*eps;
numgrad(i) = (J(theta+delta)-J(theta-delta))/(eps*);
end
m一大，eye就超出了matlab的限制了，所以，我用的上一个。
这一步句，巨慢无比，因为要进行mx2次的J函数的运算。。。。
4.结果
