UFLDL深度学习笔记 （一）反向传播与稀疏自编码

UFLDL深度学习笔记 （一）基本知识与稀疏自编码

前言

　　近来正在系统研究一下深度学习，做为新入门者，为了更好地理解、交流，准备把学习过程总结记录下来。最开始的规划是先学习理论推导；而后学习一两种开源框架；第三是进阶调优、加速技巧。越日后越要带着工做中的实际问题去作，而不能是空中楼阁式沉迷在理论资料的旧数据中。深度学习领域大牛吴恩达(Andrew Ng)老师的UFLDL教程 (Unsupervised Feature Learning and Deep Learning)提供了很好的基础理论推导，文辞既系统完备，又明白晓畅，最难得的是能把在斯坦福大学中的教学推广给全世界，尽心力而广育人，实教授之典范。
　　这里的学习笔记不是为了重复UFLDL中的推导过程，而是为了一方面补充我没有很快看明白的理论要点，另外一方面基于我对课后练习的matlab实现，记录讨论卡壳时间较长的地方。也方便其余学习者，学习过程当中UFLDL全部留白的代码模块全是本身编写的，没有查网上的代码，通过一番调试结果都达到了练习给出的参考标准，并且都按照矩阵化实现(UFLDL称为矢量化 vectorization)，代码连接见此https://github.com/codgeek/deeplearning，因此各类matlab实现细节、原因也会比较清楚，此外从实现每次练习的代码commit中能够清楚看到修改了那些代码，欢迎你们一块儿讨论共同进步。

理论推导中的要点

先上一个经典的自编码神经网络结构图。推导中经常使用的符号表征见于此

上图的神经网络结构中，输入单个向量\(\ x\\),逐层根据网络参数矩阵\(\ W, b\\)计算前向传播，最后输出向量$ h_{W,b}(x)$, 这样单个向量的偏差代价函数是

\[J{(W,b,x,y)}=\frac{1}{2}{||(h_{W,b}(x)-y||}^2\]

对全部数据代价之和即

\[\begin{align} J(W,b) &= \frac 1 m \left [\sum_{i=1}^m J{(W,b;x^{(i)},y^{(i)})} \right] +\frac \lambda 2 \sum_{l=1}^{n_l-1}\sum_{i=1}^{s_l}\sum_{j=1}^{s_{l+1}}{W_{ij}^{(l)}}^2 \\ &= \frac 1 m \left [\sum_{i=1}^m{\frac 1 2||(h_{W,b}(x^{(i)})-y^{(i)}||}^2 \right] +\frac \lambda 2 \sum_{l=1}^{n_l-1}\sum_{i=1}^{s_l}\sum_{j=1}^{s_{l+1}}{W_{ij}^{(l)}}^2 \end{align} \]

偏导与反向传导的来龙去脉

神经网络的最基础的理论是反向传导算法，反向传导的主语是偏差代价函数\(\ J_{W,b}\\)对网络系数(W，b)的偏导数。那为何要求偏导呢? 这里就引出了最优化理论，具体的问题通常被建模成许多自变量决定的代价函数，最优化的目标就是找到使代价函数最小化的自变量取值。数学方法告诉咱们令偏导为0对应的自变量取值就是最小化代价的点，不幸的是多数实际问题对用偏导为0方程求不出闭合的公式解，神经网络也不例外。为了解决这类问题，可使用优化领域的梯度降低法，梯度是对每一个自变量的偏导组成的向量或矩阵，梯度降低直白地说就是沿着对每一个自变量偏导的负方向改变自变量，降低指的是负方向，直到梯度值接近0或代价函数变化很小，即中止查找认为找到了达到最优化的全部自变量数值。 因此为了最小化稀疏自编码偏差代价函数，咱们须要求出J(W,b)对神经网络各层系数W、b的偏导数组成的梯度矩阵。

接下来在反向传导算法文章中给出了详细的推导过程，推导的结果即是残差、梯度值从输出层到隐藏层的反向传导公式

输出层的对输出值\(\ z\\)第i个份量\(\ z_i^{n_l}\\)的偏导为

\[ \delta_i^{(n_l)}= \frac {\partial} {\partial z_i^{(n_l)}} {\frac 1 2||(h_{W,b}(x^{(i)})-y^{(i)}||}^2 = -(y_i - a_i^{(n_l)}) \cdot f'( z_i^{(n_l)}) \]

后向传导时，第l+1层传导给第l层偏导值为 \[\delta_i^{(l)} = (\sum_{j=1}^{s_l+1}W_{ij}^{(l)}\delta_j^{(l+1)}) f'( z_i^{(l)}) \]

J对输出值z的偏导怎么获得对网络权值参数的偏导？

咱们知道\(\ z_i^{l+1}\\)与\(\ W,b\\)的关系\(\ z_i^{(l+1)}=\sum_{j=1}^{n_l} ({W_{i,j}^{(l)}\cdot{a_j^{(l)}+b}})\\)，运用求导链式法则可得：

\[\begin{align} \frac {\partial J_{W,b;x,y}}{\partial W_{i,j}^{(l)}} &= \frac {\partial J_{W,b;x,y}}{\partial z_i^{(l+1)}} \cdot \frac {\partial z_i^{(l+1)}}{\partial W_{i,j}^{(l)}} \\ &=\delta_i^{l+1}\cdot a_j^l\end{align}\]
原文中跳过了这一步，基于此，才能获得偏差函数对参数\(\ W, b\\)每一项的偏导数。

\[ \frac {\partial J_{W,b;x,y}}{\partial W_{i,j}^{(l)}} = \delta_i^{l+1}\cdot a_j^l\]

\[ \frac {\partial J_{W,b;x,y}}{\partial b_{i}^{(l)}} = \delta_i^{l+1}\]

结合稀疏自编码问题

前面的反向传导是神经网络的通常形式，具体到稀疏自编码，有两个关键字“稀疏”、“自”。“自”是指目标值\(\ y\\)和输入向量\(\ x\\)相等，“稀疏”是指过程要使隐藏层L1的大部分激活值份量接近0，个别份量显著大于0.因此稀疏自编码的代价函数中\(\ y\\)直接使用\(\ x\\)的值。同时加上稀疏性惩罚项,详见稀疏编码一节。
\[J_sparse(W,b)=J(W,b)+\beta\sum_{j=1}^{S_2}KL(\rho||\hat\rho_j)\]
惩罚项只施加在隐藏层上，对偏导项作对应项添加，也只有隐藏层的偏导增长一项。

\[\delta_i^{(2)} = (\sum_{j=1}^{s_2}W_{ij}^{(2)}\delta_j^{(3)}+\beta (- \frac \rho {\hat \rho_i} + \frac {1-\rho} {1-\hat \rho_i}) ) f'( z_i^{(2)}) \]

后向传导以及对网络权值\(\ W, b\\)的梯度公式不变。

在反向传导的来龙去脉段落，咱们已经知道了要用梯度降低法搜索使代价函数最小的自变量\(\ W, b\\)。实际上梯度降低也不是简单减梯度，而是有降低速率参数\(\ \alpha\\)的L-BFGS，又牵扯到一个专门的优化。为了简化要作的工做量，稀疏自编码练习中已经帮咱们集成了一个第三方实现，在minFunc文件夹中，咱们只须要提供代价函数、梯度计算函数就能够调用minFunc实现梯度降低，获得最优化参数。后续的工做也就是要只须要补上sampleIMAGES.m, sparseAutoencoderCost.m, computeNumericalGradient.m的"YOUR CODE HERE"的代码段。

exercise代码实现难点

UFLDL给你们的学习模式很到家，把周边的结构性代码都写好了matlab代码与注释，尽可能给学习者减负。系数自编码中主m文件是train.m。先用实现好的代价、梯度模块调用梯度检验，而后将上述代价、梯度函数传入梯度降低优化minFunc。知足迭代次数后退出，获得训练好的神经网络系数矩阵，补全所有待实现模块的完整代码见此，https://github.com/codgeek/deeplearning.
其中主要过程是

训练数据

训练数据来源是图片，第一步要在sampleIMAGES.m中将其转化为神经网络\(\ x\\)向量。先看一下训练数据的样子吧

sampleIMAGES.m模块是能够被后续课程复用的，尽可能写的精简通用些。从一幅图上取patchsize大小的子矩阵，要避免每一个图上取的位置相同，也要考虑并行化，循环每张图片取patch效率较低。下文给出的方法是先随机肯定行起始行列坐标，计算出numpatches个这样的坐标值，而后就能每一个patch不关联地取出来，才能运用parfor作并行循环。

sampleIMAGES.m
function patches = sampleIMAGES(patchsize, numpatches)
    load IMAGES;    % load images from disk dataSet=IMAGES;
    dataSet = IMAGES;
    patches = zeros(patchsize*patchsize, numpatches); % 初始化数组大小，为并行循环开辟空间
    [row, colum, numPic] = size(dataSet);
    rowStart = randi(row - patchsize + 1, numpatches, 1);% 从一幅图上取patchsize大小子矩阵，只须要肯定随机行起始点，范围是[1,row - patchsize + 1]
    columStart = randi(colum - patchsize + 1, numpatches, 1);% 肯定随机列起始点，范围是colum - patchsize + 1
    randIdx = randperm(numPic); % 肯定从哪一张图上取子矩阵，打乱排列，防止生成的patch顺序排列。
    parfor r=1:numpatches % 肯定了起始坐标后，每一个patch不关联，能够并行循环取
        patches(:,r) = reshape(dataSet(rowStart(r):rowStart(r) + patchsize - 1, columStart(r):columStart(r) + patchsize - 1, randIdx(floor((r-1)*numPic/numpatches)+1)),[],1);
    end
    patches = normalizeData(patches)
end

后向传播模型的矩阵实现

稀疏自编码最重要模块是计算代价、梯度矩阵的sparseAutoencoderCost.m。输入\(\ W, b\\)与训练数据data，外加稀疏性因子、稀疏项系数、权重\(\ W\\)系数，值得关注的是后向传导时隐藏层输出的残差delta2是由隐藏层与输出层的网络参数\(\ W2\\)和输出层的残差delta3决定的，不是输入层与隐藏层的网络参数\(\ W1\\)，这里最开始写错了，耽误了一些时间才调试正确。下文结合代码用注释的形式解释每一步具体做用。
sparseAutoencoderCost.m
function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ...
lambda, sparsityParam, beta, data)
%% 一维向量重组为便于矩阵计算的神经网络系数矩阵。
W1 = reshape(theta(1:hiddenSizevisibleSize), hiddenSize, visibleSize);
W2 = reshape(theta(hiddenSizevisibleSize+1:2hiddenSizevisibleSize), visibleSize, hiddenSize);
b1 = theta(2hiddenSizevisibleSize+1:2hiddenSizevisibleSize+hiddenSize);
b2 = theta(2hiddenSizevisibleSize+hiddenSize+1:end);

%% 前向传导，W1矩阵为hiddenSize×visibleSize，训练数据data为visibleSize×N_samples，训练全部数据的过程正好是矩阵相乘$\ W1*data\ $。注意全部训练数据都共用系数$\ b\ $，而单个向量的每一个份量对用使用$\ b\ $的对应份量，b1*ones(1,m)是将列向量复制m遍，组成和$\ W1*data\ $相同维度的矩阵。
    [~, m] = size(data); % visibleSize×N_samples， m=N_samples
    a2 = sigmoid(W1*data + b1*ones(1,m));% active value of hiddenlayer: hiddenSize×N_samples
    a3 = sigmoid(W2*a2 + b2*ones(1,m));% output result: visibleSize×N_samples
    diff = a3 - data; % 自编码也就意味着将激活值和原始训练数据作差值。
    penalty = mean(a2, 2); %  measure of hiddenlayer active: hiddenSize×1
    residualPenalty =  (-sparsityParam./penalty + (1 - sparsityParam)./(1 - penalty)).*beta; % penalty factor in residual error delta2
    % size(residualPenalty)
    cost = sum(sum((diff.*diff)))./(2*m) + ...
        (sum(sum(W1.*W1)) + sum(sum(W2.*W2))).*lambda./2 + ...
        beta.*sum(KLdivergence(sparsityParam, penalty));
   
    % 后向传导过程，隐藏层残差须要考虑稀疏性惩罚项，公式比较清晰。
     delta3 = -(data-a3).*(a3.*(1-a3)); %  visibleSize×N_samples
     delta2 = (W2'*delta3 + residualPenalty*ones(1, m)).*(a2.*(1-a2)); %  hiddenSize×N_samples. !!! => W2'*delta3 not W1'*delta3
     % 前面已经推导出代价函数对W2的偏导，矩阵乘法里包含了公式中l层激活值a向量与1+1层残差delta向量的点乘。
     W2grad = delta3*a2'; % ▽J(L)=delta(L+1,i)*a(l,j). sum of grade value from N_samples is got by matrix product visibleSize×N_samples * N_samples× hiddenSize. so mean value is caculated by "/N_samples"
     W1grad = delta2*data';% matrix product  visibleSize×N_samples * N_samples×hiddenSize
     
     b1grad = sum(delta2, 2);
     b2grad = sum(delta3, 2);

    % 对m个训练数据取平均
    W1grad=W1grad./m + lambda.*W1;
    W2grad=W2grad./m + lambda.*W2;
    b1grad=b1grad./m;
    b2grad=b2grad./m;% mean value across N_sample: visibleSize ×1

    % 矩阵转列向量
    grad = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)];

end

最终结果

在tarin.m中，将sparseAutoencoderCost.m传入梯度优化函数minFunc，通过迭代训练出网络参数，一般见到的各个方向的边缘检测图表示的是权重矩阵对每一个隐藏层激活值的特征，当输入值为对应特征值时，每一个激活值会有最大响应，同时的其他隐藏节点处于抑制状态。因此只须要把W1矩阵每一行的64个向量还原成8*8的图片patch，也就是特征值了，每一个隐藏层对应一个，总共100个。结果以下图。

接下来咱们验证一下隐藏层对应的输出是否知足激活、抑制的要求，将上述输入值用参数矩阵传导到隐藏层，也就是\(\ W1*W1'\\).可见，每一个输入对应的隐藏层只有一个像素是白色，表示接近1，其他是暗色调接近0，知足了稀疏性要求，并且激活的像素位置是顺序排列的。

继续改变不一样的输入层单元个数，隐藏层单元个数，能够获得更多有意思的结果！