CNN卷积神经网络

clear all; close all; clc;
addpath('../data');
addpath('../util');
load mnist_uint8;
train_x = double(reshape(train_x',28,28,60000))/255;
test_x = double(reshape(test_x',28,28,10000))/255;
train_y = double(train_y');
test_y = double(test_y');
cnn.layers = {
    struct('type', 'i') %input layer
    struct('type', 'c', 'outputmaps', 6, 'kernelsize', 5) %convolution layer
    struct('type', 's', 'scale', 2) %sub sampling layer
    struct('type', 'c', 'outputmaps', 12, 'kernelsize', 5) %convolution layer
    struct('type', 's', 'scale', 2) %subsampling layer
};
cnn = cnnsetup(cnn, train_x, train_y);
opts.alpha = 1;
opts.batchsize = 50; 
opts.numepochs = 10;
cnn = cnntrain(cnn, train_x, train_y, opts);
[er, bad] = cnntest(cnn, test_x, test_y);
plot(cnn.rL);
disp([num2str(er*100) '% error']);

上面是执行整个卷积神经网络的demo起点，主要注意cnn.layers的设置，在创建卷积神经网络时，需要根据相关的参数进行一些计算。然后是opts结构体的一些参数，alpha记录权值更新的状态，而numepochs表示进行十次迭代，battchsize表示每一批训练的数据包含50个。具体参数的作用可以参照http://www.math.duke.edu/~jvb/papers/cnn_tutorial.pdf文档的解释。

function net = cnnsetup(net, x, y)
    inputmaps = 1; 
    mapsize = size(squeeze(x(:, :, 1)));
    for l = 1 : numel(net.layers) 
        if strcmp(net.layers{l}.type, 's')
            mapsize = mapsize / net.layers{l}.scale;
            assert(all(floor(mapsize)==mapsize), ['Layer ' num2str(l) ' size must be integer. Actual: ' num2str(mapsize)]);
            for j = 1 : inputmaps
                net.layers{l}.b{j} = 0;
            end
        end
        if strcmp(net.layers{l}.type, 'c')
            mapsize = mapsize - net.layers{l}.kernelsize + 1;
            fan_out = net.layers{l}.outputmaps * net.layers{l}.kernelsize ^ 2;%本层的输出
            for j = 1 : net.layers{l}.outputmaps  %  output map
                fan_in = inputmaps * net.layers{l}.kernelsize ^ 2;
                for i = 1 : inputmaps
                    net.layers{l}.k{i}{j} = (rand(net.layers{l}.kernelsize) - 0.5) * 2 * sqrt(6 / (fan_in + fan_out));%init the weight of each map connection
                end
                net.layers{l}.b{j} = 0;
            end
            inputmaps = net.layers{l}.outputmaps;
        end
    end
    fvnum = prod(mapsize) * inputmaps;
    onum = size(y, 1);

    net.ffb = zeros(onum, 1);
    net.ffW = (rand(onum, fvnum) - 0.5) * 2 * sqrt(6 / (onum + fvnum));
end

inputmap用于计数输入map。对于输入层inputmap等于1，mapsize存放的则是输入的映射的大小，输入层的输入就是训练数据x的尺寸大小。根据层次是属于卷基层还是下采样层进行参数的计算。如果属于下采样层，则需要根据上一层的输入除以本层的scale就可以得到本层输出映射的大小，至于输出映射的个数就等于上一层卷积输入的映射数目。另外还需要根据上一层输入的映射设置偏移b等于0。对于卷积层则不相同，首先卷基层的输出映射是人为定义的，其次卷积层除掉偏移之外还需要一个卷积权重，这个权重不仅依赖与输出映射，还依赖于输入映射。首先计算的是mapsize的大小，根据输出层大小以及卷积核的大小和卷积平移的步长来计算的，这里平移步长为1。首先计算fan_out和fan_in用作归一化处理。下面的双重循环就是计算权重，每一个输入映射都对应outputmaps个输出映射，至于每一个输出组的偏移是可以共享的。最后计算的是输出层的权重和偏移。最后输出层将最后得到的输出映射以及映射大小转化为一个向量，向量的大小为fvnum。然后根据输入的标签得到最终所需要的标签特征的向量大小为onum。因此最后的偏移是onum个，而w的权值是onum*fvnum。因为卷积神经网络本质上是一个分类器，所以onum实际代表的是类标签的数目，而fvnum则表示根据卷积神经网络得到的一个x的特征，经过分类处理将x判断为onum个标签中的某一个。

function net = cnntrain(net, x, y, opts)
    m = size(x, 3);
    numbatches = m / opts.batchsize;
    if rem(numbatches, 1) ~= 0
        error('numbatches not integer');
    end
    net.rL = [];
    for i = 1 : opts.numepochs
        disp(['epoch ' num2str(i) '/' num2str(opts.numepochs)]);
        tic;
        kk = randperm(m);
        for l = 1 : numbatches
            batch_x = x(:, :, kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));
            batch_y = y(:,    kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));
            net = cnnff(net, batch_x);
            net = cnnbp(net, batch_y);
            net = cnnapplygrads(net, opts);
            if isempty(net.rL)
                net.rL(1) = net.L;
            end
            net.rL(end + 1) = 0.99 * net.rL(end) + 0.01 * net.L;
        end
        toc;
    end
end

总体的训练过程没什么可说的，主要是根据传递进来的数据随机生成批次训练的数据，然后调用下一层次的前向传播和反向传播等等操作。最后 根据新得到的损失值来更新原来的损失值。

function net = cnnff(net, x)
    n = numel(net.layers);
    net.layers{1}.a{1} = x;
    inputmaps = 1;
    for l = 2 : n%  for each layer
        if strcmp(net.layers{l}.type, 'c')
            for j = 1 : net.layers{l}.outputmaps   %  for each output map
                z = zeros(size(net.layers{l - 1}.a{1}) - [net.layers{l}.kernelsize - 1 net.layers{l}.kernelsize - 1 0]);
                for i = 1 : inputmaps   %  for each input map
                    z = z + convn(net.layers{l - 1}.a{i}, net.layers{l}.k{i}{j}, 'valid');
                end
                net.layers{l}.a{j} = sigm(z + net.layers{l}.b{j});
            end
            inputmaps = net.layers{l}.outputmaps;
        elseif strcmp(net.layers{l}.type, 's')
            for j = 1 : inputmaps
                z = convn(net.layers{l - 1}.a{j}, ones(net.layers{l}.scale) / (net.layers{l}.scale ^ 2), 'valid');   %  !! replace with variable
                net.layers{l}.a{j} = z(1 : net.layers{l}.scale : end, 1 : net.layers{l}.scale : end, :);
            end
        end
    end
    net.fv = [];
    for j = 1 : numel(net.layers{n}.a)
        sa = size(net.layers{n}.a{j});
        net.fv = [net.fv; reshape(net.layers{n}.a{j}, sa(1) * sa(2), sa(3))];
    end
    net.o = sigm(net.ffW * net.fv + repmat(net.ffb, 1, size(net.fv, 2)));
end

首先函数根据传递进来的参数得到层次的数目，然后将输入的训练样本整个传递给输入层，设置初始的输入映射为1，然后根据所属的层训练样本。在卷积层中，根据上一层的输出得到本层的输入，然后根据本层卷积数目以及步长计算本层的输出的大小Z，并将z初始化为0.这里需要注意的是z为三维矩阵，所以在下面的convn卷积函数中可以一次性处理一批50个数据，并且由于参数为valid，所以得到的矩阵为size（a）-size（k）+1。最后根据激活函数计算本层的输出。如果是在下采样层，则直接进行卷积，然后根据结果进行下采样，不过这里的convn最后一个参数应该是same。输出层则直接将参数转化为以为数组，最后sa（3）实际上是一批训练样本的数目。最后根据激活函数计算输出的类别，net.o的大小为（onum*sa（3））。

function net = cnnbp(net, y)
    n = numel(net.layers);
    net.e = net.o - y;
    net.L = 1/2* sum(net.e(:) .^ 2) / size(net.e, 2);
    net.od = net.e .* (net.o .* (1 - net.o));   %  output delta
    net.fvd = (net.ffW' * net.od);              %  feature vector delta
    if strcmp(net.layers{n}.type, 'c')         %  only conv layers has sigm function
        net.fvd = net.fvd .* (net.fv .* (1 - net.fv));
    end
    sa = size(net.layers{n}.a{1});
    fvnum = sa(1) * sa(2);
    for j = 1 : numel(net.layers{n}.a)
        net.layers{n}.d{j} = reshape(net.fvd(((j - 1) * fvnum + 1) : j * fvnum, :), sa(1), sa(2), sa(3));
    end
    for l = (n - 1) : -1 : 1
        if strcmp(net.layers{l}.type, 'c')
            for j = 1 : numel(net.layers{l}.a)
                net.layers{l}.d{j} = net.layers{l}.a{j} .* (1 - net.layers{l}.a{j}) .* (expand(net.layers{l + 1}.d{j}, [net.layers{l + 1}.scale net.layers{l + 1}.scale 1]) / net.layers{l + 1}.scale ^ 2);
            end
        elseif strcmp(net.layers{l}.type, 's')
            for i = 1 : numel(net.layers{l}.a)
                z = zeros(size(net.layers{l}.a{1}));
                for j = 1 : numel(net.layers{l + 1}.a)
                     z = z + convn(net.layers{l + 1}.d{j}, rot180(net.layers{l + 1}.k{i}{j}), 'full');
                end
                net.layers{l}.d{i} = z;
            end
        end
    end
    %%  calc gradients
    for l = 2 : n
        if strcmp(net.layers{l}.type, 'c')
            for j = 1 : numel(net.layers{l}.a)
                for i = 1 : numel(net.layers{l - 1}.a)
                    net.layers{l}.dk{i}{j} = convn(flipall(net.layers{l - 1}.a{i}), net.layers{l}.d{j}, 'valid') / size(net.layers{l}.d{j}, 3);
                end
                net.layers{l}.db{j} = sum(net.layers{l}.d{j}(:)) / size(net.layers{l}.d{j}, 3);
            end
        end
    end
    net.dffW = net.od * (net.fv)' / size(net.od, 2);
    net.dffb = mean(net.od, 2);
    function X = rot180(X)
        X = flipdim(flipdim(X, 1), 2);
    end
end

根据已有的参数计算出误差net.e，并且统计得到最小的误差——也就是损失函数net.L。那么根据求导法则不难求出net.od，其中维1/2(X^2)外加sigmoid函数的求导公式就可以推导出来。注意这里是对偏移b求导，所以得到的是传播误差，也就是note中所说的delta。然后根据note上面的传播公式，可以得到输出层的上一层的传播函数，当上一层是卷积层则需要进行传播。传播公式为note里面的公式4。公式4的数学推导过程如下：

delta         E/u

-------- =  --------

delta1      E/u1

根据上述公式进行数学推导就可以得到note中的公式4.

在第一个循环中将输出层的灵敏度也就是delta进行分解，并往上层传递得到输出层上一层的输出的delta。假设第l-1层是卷积层，则l层为下采样层，所以对delta的传递需要分为两次处理：首先，将下层的delta进行一个上采样，这个在expand函数里面执行；第二，根据delta传递公式往上层传递delta值，其中对本层的求导公式就是最前面的那一部分。如果l-1层是下采样曾，则其l层为卷积层。正好卷积层的输入等于下采样层的输出；所以传递过程仅仅是卷积层的互相关操作就可以了。由于自相关是一个卷积过程，所以需要有一个累加。到这里BP操作基本完成，从整个过程中可以看到，误差会从输出层一直传递到上层的输入层，这就是深度神经网络相对于原来的神经网络的优势，误差不会在BP过程中丢失。

在计算得到误差delta之后需要更新参数，由于在下采样层没有参数，所以下采样曾不需要计算参数，而仅仅只需要计算卷积层的参数。根据计算梯度的公式可以得到卷积核和偏移的求导方式。最后更新输出层权值和偏移的数值。

function net = cnnapplygrads(net, opts)
    for l = 2 : numel(net.layers)
        if strcmp(net.layers{l}.type, 'c')
            for j = 1 : numel(net.layers{l}.a)
                for ii = 1 : numel(net.layers{l - 1}.a)
                    net.layers{l}.k{ii}{j} = net.layers{l}.k{ii}{j} - opts.alpha * net.layers{l}.dk{ii}{j};
                end
                net.layers{l}.b{j} = net.layers{l}.b{j} - opts.alpha * net.layers{l}.db{j};
            end
        end
    end
    net.ffW = net.ffW - opts.alpha * net.dffW;
    net.ffb = net.ffb - opts.alpha * net.dffb;
end

最后一个函数就是权值的更新，通过上面的分析，没什么难度。上面就是整个CNN卷积神经网络的训练过程。