Add sparse CNN

SCNN/scnnapplygrads.m

@@ -0,0 +1,15 @@
+function net = scnnapplygrads(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

SCNN/scnnbp.m

@@ -0,0 +1,64 @@
+function net = scnnbp(net, y)
+n = numel(net.layers);
+
+%   error
+net.e = net.o - y;
+%  loss function
+net.L = 1/2* sum(net.e(:) .^ 2) / size(net.e, 2);
+
+%%  backprop deltas
+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
+
+%  reshape feature vector deltas into output map style
+sa = sptsize(net.layers{n}.a{1});
+fvnum = sa(1) * sa(2);
+for j = 1 : numel(net.layers{n}.a)
+    net.layers{n}.d{j} = tensor2spcell(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)
+            expand_param =  [net.layers{l + 1}.scale net.layers{l + 1}.scale];
+            expanded_delta = spcell_expand(net.layers{l + 1}.d{j},expand_param);
+            net.layers{l}.d{j} = spcell_mul(spcell_mul(net.layers{l}.a{j},spcell_sub(1,net.layers{l}.a{j})),...
+                spcell_div(expanded_delta,net.layers{l + 1}.scale ^ 2));
+        end
+    elseif strcmp(net.layers{l}.type, 's')
+        for i = 1 : numel(net.layers{l}.a)
+            z = sptzeros(sptsize(net.layers{l}.a{1}));
+            for j = 1 : numel(net.layers{l + 1}.a)
+                z = spcell_add(z,spconvn(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)
+            n_iter_i = numel(net.layers{l - 1}.a);
+            spcell_temp = cell(1,n_iter_i);
+            parfor i = 1 : n_iter_i
+                spcell_temp{i} =  spconvn(sptflipall(net.layers{l - 1}.a{i}),net.layers{l}.d{j}, 'valid') / size(net.layers{l}.d{j}, 3);
+            end
+            for i = 1:n_iter_i
+                net.layers{l}.dk{i}{j} = spcell_temp{i};
+            end
+            net.layers{l}.db{j} = sptsum(net.layers{l}.d{j}) / numel(net.layers{l}.d{j});
+        end
+    end
+end
+net.dffW = net.od * (net.fv)' / size(net.od, 2);
+net.dffb = mean(net.od, 2);
+
+    function X = rot180(X)
+        X = flip(flip(X, 1), 2);
+    end
+end

SCNN/scnnff.m

@@ -0,0 +1,46 @@
+function net = scnnff(net, x)
+%% DOCUMENT
+%{
+type of x: cell_array<sparse<double>(x,y)>(n)
+%}
+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')
+        %  !!below can probably be handled by insane matrix operations
+        for j = 1 : net.layers{l}.outputmaps   %  for each output map
+            %  create temp output map
+            z = sptzeros(sptsize(net.layers{l - 1}.a{1}) - [net.layers{l}.kernelsize - 1 net.layers{l}.kernelsize - 1 0]);
+            temp_matrix = cell(1,inputmaps);
+            for i = 1 : inputmaps   %  for each input map
+                %  convolve with corresponding kernel and add to temp output map
+                temp_matrix{i} = spconvn(net.layers{l - 1}.a{i}, net.layers{l}.k{i}{j}, 'valid');
+            end
+            for i = 1 : inputmaps 
+                z = spcell_add(z,temp_matrix{i});
+            end
+            %  add bias, pass through nonlinearity
+            net.layers{l}.a{j} = spcell_sigm(spcell_add(z,net.layers{l}.b{j}));
+        end
+        %  set number of input maps to this layers number of outputmaps
+        inputmaps = net.layers{l}.outputmaps;
+    elseif strcmp(net.layers{l}.type, 's')
+        %  downsample
+        for j = 1 : inputmaps
+            z = spconvn(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} = sptsubsample(z,net.layers{l}.scale);
+        end
+    end
+end
+
+%  concatenate all end layer feature maps into vector
+net.fv = [];
+for j = 1 : numel(net.layers{n}.a)
+    % sa = sptsize(net.layers{n}.a{j});
+    net.fv = [net.fv; sptflatten(net.layers{n}.a{j})];
+end
+%  feedforward into output perceptrons
+net.o = sigm(net.ffW * net.fv + repmat(net.ffb, 1, size(net.fv, 2)));
+end

SCNN/scnnnumgradcheck.m

@@ -0,0 +1,79 @@
+function scnnnumgradcheck(net, x, y)
+    epsilon = 1e-4;
+    er      = 1e-8;
+    n = numel(net.layers);
+    for j = 1 : numel(net.ffb)
+        net_m = net; net_p = net;
+        net_p.ffb(j) = net_m.ffb(j) + epsilon;
+        net_m.ffb(j) = net_m.ffb(j) - epsilon;
+        net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
+        net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
+        d = (net_p.L - net_m.L) / (2 * epsilon);
+        e = abs(d - net.dffb(j));
+        if e > er
+            error('numerical gradient checking failed');
+        end
+    end
+
+    for i = 1 : size(net.ffW, 1)
+        for u = 1 : size(net.ffW, 2)
+            net_m = net; net_p = net;
+            net_p.ffW(i, u) = net_m.ffW(i, u) + epsilon;
+            net_m.ffW(i, u) = net_m.ffW(i, u) - epsilon;
+            net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
+            net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
+            d = (net_p.L - net_m.L) / (2 * epsilon);
+            e = abs(d - net.dffW(i, u));
+            if e > er
+                error('numerical gradient checking failed');
+            end
+        end
+    end
+
+    for l = n : -1 : 2
+        if strcmp(net.layers{l}.type, 'c')
+            for j = 1 : numel(net.layers{l}.a)
+                net_m = net; net_p = net;
+                net_p.layers{l}.b{j} = net_m.layers{l}.b{j} + epsilon;
+                net_m.layers{l}.b{j} = net_m.layers{l}.b{j} - epsilon;
+                net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
+                net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
+                d = (net_p.L - net_m.L) / (2 * epsilon);
+                e = abs(d - net.layers{l}.db{j});
+                if e > er
+                    error('numerical gradient checking failed');
+                end
+                for i = 1 : numel(net.layers{l - 1}.a)
+                    for u = 1 : size(net.layers{l}.k{i}{j}, 1)
+                        for v = 1 : size(net.layers{l}.k{i}{j}, 2)
+                            net_m = net; net_p = net;
+                            net_p.layers{l}.k{i}{j}(u, v) = net_p.layers{l}.k{i}{j}(u, v) + epsilon;
+                            net_m.layers{l}.k{i}{j}(u, v) = net_m.layers{l}.k{i}{j}(u, v) - epsilon;
+                            net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
+                            net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
+                            d = (net_p.L - net_m.L) / (2 * epsilon);
+                            e = abs(d - net.layers{l}.dk{i}{j}(u, v));
+                            if e > er
+                                error('numerical gradient checking failed');
+                            end
+                        end
+                    end
+                end
+            end
+        elseif strcmp(net.layers{l}.type, 's')
+%            for j = 1 : numel(net.layers{l}.a)
+%                net_m = net; net_p = net;
+%                net_p.layers{l}.b{j} = net_m.layers{l}.b{j} + epsilon;
+%                net_m.layers{l}.b{j} = net_m.layers{l}.b{j} - epsilon;
+%                net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
+%                net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
+%                d = (net_p.L - net_m.L) / (2 * epsilon);
+%                e = abs(d - net.layers{l}.db{j});
+%                if e > er
+%                    error('numerical gradient checking failed');
+%                end
+%            end
+        end
+    end
+%    keyboard
+end

SCNN/scnnsetup.m

@@ -0,0 +1,36 @@
+function net = scnnsetup(net, x, y)
+    %assert(~isOctave() || compare_versions(OCTAVE_VERSION, '3.8.0', '>='), ['Octave 3.8.0 or greater is required for CNNs as there is a bug in convolution in previous versions. See http://savannah.gnu.org/bugs/?39314. Your version is ' myOctaveVersion]);
+    inputmaps = 1;
+    mapsize = size(x{1});
+
+    for l = 1 : numel(net.layers)   %  layer
+        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  %  input map
+                    net.layers{l}.k{i}{j} = (rand(net.layers{l}.kernelsize) - 0.5) * 2 * sqrt(6 / (fan_in + fan_out));
+                end
+                net.layers{l}.b{j} = 0;
+            end
+            inputmaps = net.layers{l}.outputmaps;
+        end
+    end
+    % 'onum' is the number of labels, that's why it is calculated using size(y, 1). If you have 20 labels so the output of the network will be 20 neurons.
+    % 'fvnum' is the number of output neurons at the last layer, the layer just before the output layer.
+    % 'ffb' is the biases of the output neurons.
+    % 'ffW' is the weights between the last layer and the output neurons. Note that the last layer is fully connected to the output layer, that's why the size of the weights is (onum * fvnum)
+    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

SCNN/scnntest.m

@@ -0,0 +1,9 @@
+function [er, bad] = scnntest(net, x, y)
+    %  feedforward
+    net = scnnff(net, x);
+    [~, h] = max(net.o);
+    [~, a] = max(y);
+    bad = find(h ~= a);
+
+    er = numel(bad) / size(y, 2);
+end

SCNN/scnntrain.m

@@ -0,0 +1,35 @@
+function net = scnntrain(net, x, y, opts)
+    m = numel(x);
+    numbatches = m / opts.batchsize;
+    if rem(numbatches, 1) ~= 0
+        error('numbatches not integer');
+    end
+    net.rL = [];
+    wb = waitbar(0,'Training...(Sparse)');
+    sum_count = opts.numepochs * numbatches;
+    t_index = 1;
+    for i = 1 : opts.numepochs
+        disp(['epoch ' num2str(i) '/' num2str(opts.numepochs)]);
+        tic;
+        kk = randperm(m);
+        for l = 1 : numbatches
+
+            waitbar(t_index/sum_count,wb,...
+                sprintf('Training...Epoch:%d of %d  BatchIndex%d of %d',i,opts.numepochs,l,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 = scnnff(net, batch_x);
+            net = scnnbp(net, batch_y);
+            net = scnnapplygrads(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;
+            t_index = t_index + 1;
+        end
+        toc;
+    end
+
+end

SCNN/spcell_add.m

@@ -0,0 +1,17 @@
+function x = spcell_add(x, y)
+%{
+x,y,z type: array<sparse<double>(x,y)>(n)
+%}
+if iscell(y) && iscell(x)
+    n = numel(x);
+    for i = 1:n
+        x{i} = x{i} + y{i};
+    end
+elseif isa(y,'double') && iscell(x)
+    x = spcell_ewop(x,@(u) u+y);
+elseif isa(x,'double') && iscell(y)
+    x = spcell_ewop(y,@(u) u+x);
+else
+    x = x + y;
+end
+end

SCNN/spcell_div.m

@@ -0,0 +1,25 @@
+function x = spcell_div(x, y)
+%{
+    �������ܣ���spcell�ģ��㣩��������
+%}
+    if iscell(y) && iscell(x)
+        n = numel(x);
+        for i = 1:n
+            x{i} = x{i} ./ y{i};
+        end
+    elseif isa(y,'double') && iscell(x)
+        % x = spcell_ewop(x,@(u) u./y);
+        n = numel(x);
+        for i = 1:n
+            x{i} = x{i}./y;
+        end
+    elseif isa(x,'double') && iscell(y)
+        n = numel(y);
+        for i = 1:n
+            y{i} = x./y{i};
+        end
+        x = y;
+    else
+        x = x./y;
+    end
+end

SCNN/spcell_ewop.m

@@ -0,0 +1,9 @@
+function x = spcell_ewop(x, operation)
+%{
+    �������ܣ�ʹ�ò�������ÿһ��Ԫ�ؽ��в���
+%}
+    n = numel(x);
+    for i = 1:n
+        x{i} = operation(x{i});
+    end
+end

SCNN/spcell_expand.m

@@ -0,0 +1,7 @@
+function x = spcell_expand(x, expand_param)
+    % x = spcell_ewop(x,@(u) expand(u,expand_param));
+    n = numel(x);
+    for i = 1:n
+        x{i} = expand(x{i},expand_param);
+    end
+end

SCNN/spcell_mul.m

@@ -0,0 +1,24 @@
+function x = spcell_mul(x, y)
+%{
+x,y,z type: array<sparse<double>(x,y)>(n)
+%}
+if iscell(y) && iscell(x)
+    n = numel(x);
+    for i = 1:n
+        x{i} = x{i} .* y{i};
+    end
+elseif isa(y,'double') && iscell(x)
+    n = numel(x);
+    for i = 1:n
+        x{i} = x{i}.*y;
+    end
+elseif isa(x,'double') && iscell(y)
+    n = numel(y);
+    for i = 1:n
+        y{i} = y{i}.*x;
+    end
+    x = y;
+else
+    x = x .* y;
+end
+end

SCNN/spcell_sigm.m

@@ -0,0 +1,7 @@
+function X = spcell_sigm(X)
+    % Sigmoid 激活函数
+    n = numel(X);
+    for i = 1:n
+        X{i} = 1./(1+exp(-X{i}));
+    end
+end