从软件工程的角度写机器学习6——深度学习之卷积神经网络（CNN）实现

卷积神经网络（CNN）实现

背景

卷积神经网络广泛用于图像检测，它的实现原理与传统神经网络基本上是一样的，因此将普遍意义的神经网络和卷积神经网络的实现合成一篇。

神经网络实现思路

“扔掉神经元”

尽管所有教程在介绍神经网络时都会把一大堆神经元画出来，并且以输入节点——神经元——输出结点连线，但是，在编程实现时，基于神经元去编程是低效的。典型如这篇经典文章里面的代码：
http://blog.csdn.net/zzwu/article/details/575125。
比较合适的方法是将神经网络每个层仅仅视为一个矩阵算符，对输入作变换后传递给下一层。基于矩阵运算的编程，思路清晰、容易校验，最重要的是便于后续性能优化，足够快。
因此，在写神经网络算法时，建议把“神经元”这一概念扔掉，在推导出矩阵变换公式之后，这一概念对我们工程师而言已经没有意义，我们面对的，仅仅是一个个的矩阵算符，理解算符并实现就可以了。实现神经网络，就是实现各类矩阵算符，并按顺序连接起来。

网络结构的表示

如图所示，一个神经网络Net由若干个Layer和一个全局参数矩阵Parameters（参数矩阵高为1,实则为一个向量）构成，每个Layer拥有自己独立的算符Op和运算缓存Cache，并将全局参数矩阵中的一部分映射为自己的参数矩阵P。

Layer结构

每个层由算符、参数和缓存构成。
算符负责实现矩阵变换：

Y=f(X,P)

上图是一个两层的神经网络向量变换过程。
batch size 表示一次进行计算的向量个数，input width 为输入向量维度，output width 为神经网络的输出向量维度。
算符对矩阵中的每个向量进行操作，对应地转换为另一个向量。算符实现的是向量变换的功能，之所以要用矩阵的形式表示，一方面，在随机批量梯度下降算法中，需要一次性抽取一批样本作训练，这样本身就形成矩阵。另一方面，要加大运算量，便于工程上后续作多线程/异构计算优化。多线程/异构计算的启动是有额外开销的（任务调度、kernel编译、内存传输等等），单次运算量太小会使得优化得不偿失。

Cache为缓存，仅仅做预测时，这是不需要的，但在训练过程（BP算法）中，往往需要缓存该层的输入输出，以便后续计算梯度。

Layer中的参数矩阵由网络中的全局参数矩阵截取映射而来。
对每一层，设X为输入矩阵，Y为输出矩阵，P为该层参数矩阵，则有：

Y=f(X,P)

Layer算符实现f(X,P)，Layer维护相应的cache和paramters

预测过程

预测就是一次前向传播，每一个Layer算出Y值后，作为下一层的X值传入。
设有3个Layer，那么输出结果的表示就是：

Y=f3(f2(f1(X,P1),P2),P3)

训练过程

神经网络算法是一系列矩阵算符的叠加，训练神经网络就是求出最佳参数矩阵。
这个训练过程一般基于随机梯度下降，计算梯度时采用反向传播(backward)方式。

随机梯度下降

随机梯度下降（严格来说是随机批量梯度下降）的算法描述如下：
1、从样本集中随机抽取n个样本。
2、计算这批样本对参数P所产生的梯度ΔP。
3、更新参数：P=(1−λ)P−αΔP。
4、回到第1步，循环执行iteration次。

在执行随机批量梯度下降算法时，需要设定如下超参数：
1、梯度下降的步长α
2、每次训练抽取的样本数n，也就是batch size
3、正则惩罚项λ，
4、迭代次数iteration

有些文献中，这些超参数并不是固定的，而是随着迭代次数或误差总值做变化，此处暂不考虑。

后向传播算法

设Y⎯⎯⎯为目标输出矩阵，则损失函数被定义为：

L=12||Yi−Yi⎯⎯⎯||2+12λ||P||2

λ为前面所说的正则项，在梯度下降算法中统一考虑。
经过不严格的推导，可得：

∂L∂X=∂Y∂X(Y−Y⎯⎯⎯)

∂L∂P=∂Y∂P(Y−Y⎯⎯⎯)

∂L∂P就是该层的参数梯度，求出之后先缓存，在上级的梯度下降算法中统一更新参数。
∂L∂X就是X−X⎯⎯⎯，即上一层的输出残差。
每一层求出这两个矩阵，并把∂L∂X作为上一层的输出残差Y−Y⎯⎯⎯传回去，在上一层继续求梯度，这就是后向传播算法。

输出层残差的计算

在后向传播算法中，有了最后一层的输出残差，就能逐步往前更新各层的参数，计算残差只需要将预测矩阵和目标矩阵作减法就可以。因此这个问题等同于怎么得到目标矩阵Y⎯⎯⎯。
对于回归问题，Y⎯⎯⎯中每行是一个1维向量，就是标注的一个实数值。对于自动编码器，Y⎯⎯⎯就是第一层的输入矩阵X。

对于分类问题，用Softmax为最后一层时，Y⎯⎯⎯是一个分布矩阵，每一行在标注的那一个位标1，其他元素为0。
如下图示例：

主要算符实现

前面讲述了一个通用的神经网络结构设计，现在需要到具体到每个层的实现。

卷积层（Convolution）

这个是卷积神经网络的核心，也是最难理解的一层。
英文教程参考：
http://cs231n.github.io/convolutional-networks/

卷积层、池化层都是以三维数组的方式处理矩阵中的一行，总体来说，将输入矩阵看成四维数组处理，其得到的也将是四维数组。
这是因为，CNN一般处理的是图像，图像数据原本就是3维的（宽、高、通道数），在映射为矩阵时才变为矩阵中的一行，按图像真实性质将输入数据重构为3维，可以取得良好效果。

如图所示：
输入矩阵 X 被表示为 batch size 个iw*ih*kd的立方体，batch size 为输入样本数。
参数矩阵 P 有 filter number （后面简写为kn）行，每一行是一个滤波器，它包含kh*kw*kd个系数及一个常数项C。

Y=filter(X,P)

每一个滤波器均与输入向量作一次滤波，得到一个 oh*ow 的平面，由于有kn个滤波器，得到的就是 oh*ow*kn 的输出向量。
oh和ow的计算公式中，p为输入矩阵补0的大小，s为产生输出的间隔，目前简单起见就设p=0,s=1。

滤波运算产生平面的公式如下：
设In为输入的三维数组，Out为其中一个输出平面，Kp为当前所取的滤波器，那么：

Out(oi,oj)=C+∑i=0kw∑j=0kh∑k=0kdKp(i,j,k)⋅In(oi+i,oj+j,k)

卷积层终究只是一个线性变换。计算其梯度的原则就是对该分量找到所有与它相关的参数，求和叠加。

仅考虑s=1和p=0的情况，
求输入残差ΔX，那么对X(x,y)，先将x转化为三维坐标：i,j,k，然后其值就是

ΔX(i,j,k,y)=∑p=0kn∑u=0kw∑v=0khKp(u,v,k)⋅ΔY(i−u,j−v,p,y)

注：ΔY(i,j,p,y)在i<0或j<0时取0

对于ΔP，其公式为：

ΔPp(i,j,k)=∑y=0n∑u=0ow∑v=0ohΔY(u,v,p,y)⋅X(ow+u,oh+v,k,y)

由于卷积层的运算非常大，且运算特殊，完全基于矩阵的四则运算虽能实现（如caffe的GEMM方法）但性能不是最优，建议独立为其设立矩阵算符。

池化层（Pooling）

这一层依然把输入矩阵中的一行当三维数组处理，将平面缩小，深度不变：

iw∗ih∗d−→−−Pooliws∗ihs∗d

s为缩小倍率。

计算公式可表示为

Y(i,j,k,y)=Pools,su,v=0,0X(i∗s+u,j∗s+v,k,y)

Pool为Max或Mean

池化层没有参数，只需要求输入残差。
均值池化是一个线性变换，最大池化是一个分段线性变换。
均值法的输入残差计算如下式：

ΔX(i,j,k,y)=1s2ΔY(i/s,j/s,k,y)

最大值法的输入残差计算：

ΔX(i,j,k,y)=(X(i,j,k,y)=max)?ΔY(i/s,j/s,k,y):0

内积层（InnerProduct/FullConnect）

这一层又称全连接层。因为输入向量中的每一维和输出向量中的每一维都有一个权值，因此参数个数相当多。

Y=XP

计算来看，内积层/全连接层就是一个矩阵的线性变换，其后向传播公式可以简单推得。

ΔX=ΔYP,ΔP=XTΔY

此处没有考虑常数项，考虑常数项的话把输入矩阵后面补一列1就可以了。

正则层（Relu）

这一层作用是把所有数校正为非零的。

Y=X>0?X:0

这一层没有参数，只需要计算输入残差，公式如下：

ΔX=X>0?ΔY:0

逻辑回归层（SoftMax）

公式参考：
http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92

此处设输入矩阵的宽为w
考虑到前面可以接内积层，这一层就不需要设参数了，直接做变换即可：

Y(x,y)=e−X(x,y)∑wi=0e−X(i,y)

梯度推导
此处只需要计算输入残差，经过求导之后，得到下面式子：

ΔX(x,y)=Y(x,y)(1−Y(x,y))ΔY(x,y)

简单些的表示是对矩阵中每个元素均有：

Δx=y(1−y)Δy

代码实现

Layer

算符

由于代码中打不出Δ这种符号，上面推演公式中的ΔX对应before_diff，X对应before，ΔY对应after_diff，Y对应after，P对应parameters，ΔP对应parameters_diff。

class ILayerOperator{public:    /*根据输入矩阵的宽（输入向量维度），计算本算符的输出矩阵宽(输出向量维度)*/    virtual size_t vComputeOutputWidth(size_t w) const;    /*前向传播，计算输出矩阵*/    virtual void vForward(const Matrix* before, Matrix* after/*Output*/, const Matrix* parameters) const = 0;    /*后向传播，计算输入残差和参数梯度*/    virtual void vBackward(const Matrix* after_diff, const Matrix* after, const Matrix* before, Matrix* before_diff/*Output*/, const Matrix* parameters, Matrix* parameters_diff/*Output*/) const = 0;    /*对该层所需参数的初始化算法*/    virtual size_t vInitParameters(Matrix* parameters) const = 0;    virtual ~ ILayerOperator(){}protected:    ILayerOperator(){}};

具体各Layer算符这里不再讲述。

训练用Layer

class TrainLayer{public:    //参数映射，返回映射后的偏移值    size_t mapParameters(Matrix* parameters, size_t offset);    //参数梯度目标值映射，parameters和parameters_diff同大小    size_t mapParametersDiff(Matrix* parameters_diff, size_t offset);    //前向传播，得到预测结果    Matrix* forward(Matrix* input);    //后向传播，计算本层的参数梯度和输入梯度，并将输入梯度传到上一层    double backward(Matrix* output_diff);private:    TrainLayer* mBefore;    TrainLayer* mNext;    /*在forward时，保存本层的输入输出，以便backward时使用*/    Matrix* mInputCache;    Matrix* mOutputCache;    /*参数矩阵和参数梯度矩阵的引用*/    Matrix* mParameterRef;    Matrix* mParameterDiffRef;};

预测用Layer

class PredictLayer{public:    size_t mapParameters(Matrix* parameters, size_t offset);    Matrix* forward(Matrix* input);private:    PredictLayer* mNext;//预测时只需要知道下一层    Matrix* mParameterRef;//参数引用};

训练相关

训练器

class NNLearner : public ILearner{public:    /*这里用Node表示各个层的信息，一般而言，可以写成json，然后解析json而得，在构造函数中确定默认输入向量大小，创建所有Layer的算符*/    NNLearner(Node* info);    virtual ~NNLearner();    /*这个函数所做的事情如下：    1、基于X的宽，创建各个算符的输入输出缓存，初始化参数配置，从而创建逐层相连TrainLayer，进而创建梯度计算的类NNDerivativeFunction。    2、将Y展开为目标向量，与X合并成为梯度下降所需要的混合矩阵    3、根据各个算符所需要参数的总大小，创建一个总参数矩阵，映射给TrainLayer，并用算符对其进行初始化。    4、创建一个梯度下降算法类，调节参数矩阵的值    5、最后按算符重建一系列的TestLayer，并映射参数矩阵的值，将第一个TestLayer和参数矩阵打包，即为预测器*/    virtual IPredictor* vLearn(const Matrix* X, const Matrix* Y) const;private:    /*依次存储各个layer的算符*/    std::vector<ILayerOperator*> mLayerOps;    size_t mDefaultInputWidth;};

梯度算符

class NNDerivativeFunction : public IGradientDecent::DerivativeFunction{public:    /*M为混合矩阵，对矩阵的每一行，前mOutputSize为输出向量，后面的是输入向量，在计算时先将输入矩阵X抽出来，输入mFirst前向传播，得到输出矩阵Y，然后抽出输出矩阵YP，计算残差，从mLast开始反向传播，计算完成后，输出参数残差parameters_diff*/    virtual Matrix* vCompute(Matrix* coefficient, Matrix* M) const;private:    TrainLayer* mFirst;    TrainLayer* mLast;    size_t mOutputSize;};

随机梯度下降算法

class StochasticGradientDecent : public IGradientDecent{public:    virtual void vOptimize(Matrix* coefficient, Matrix* X, const DerivativeFunction* delta, double alpha, int iteration) const    {        for (int i=0; i<iteration; ++i)        {            Matrix* selectX = Matrix::randomeSelect(X, mBatchSize);            Matrix* deltaC = delta->vCompute(coefficient, selectX);            /*更新参数： C = (1-lambda)*C-alpha*deltaC*/            Matrix::linear(coefficient, coefficient, 1.0-mLambda, deltaC.get(), -alpha);            delete deltaC;            delete selectX;        }    }private:    int mBatchSize;    double mLambda;};

预测器

class NNPredictor : public IPredictor{public:    /*Forward就可以了*/    virtual Matrix* vPredict(Matrix* X) const;private:    TestLayer* mFirst;    Matrix* mParameters;};

代码结构图如下：