CNN3 im2col

原理参照Caffe学习 五 conv_layer与im2col。

以下是一个比较，使用自己写的conv2d的卷积速度非常慢。

im2col转换进行的卷积与sg.convolve2d（不旋转180°）速度接近。

除此之外，im2col还使得反向传播更为简单。

im2col将卷积操作变为矩阵乘法，反向传播便和全连接层一样。

import scipy.signal.signaltools as sg
import numpy as np
import timeit

b=range(20*4*5*5)
b=np.reshape(b,(20,4,5,5))*0.001
a=range(10*4*28*28)
a=np.reshape(a,(10,4,28,28))*0.001


def imcol(a, v):
    n, m, h, w = np.shape(a)
    vn, vm, vh, vw = np.shape(v)
    t = h - vh + 1
    kh = t ** 2  
    kw = vh * vh  
    z = np.zeros((n, kh, kw * m))
    for i in range(n):
        for j in range(m):
            kt = []
            for it in range(t):
                for jt in range(t):
                    kt.append((a[i, j, it:it + vh, jt:jt + vh]).flatten())
            z[i, :,j*kw:j*kw+kw] = kt
    z=np.reshape(z,(n*kh,kw*m))
    v=np.reshape(v,(-1,kw*m))
    res = np.dot(z, v.T)
    res = np.reshape(res, (n, kh, -1))
    res = np.rollaxis(res, 2, 1)
    return np.reshape(res, (n, -1, t, t))

def conv2d(a, v):
    ah, aw = np.shape(a)
    vh, vw = np.shape(v)
    k = [[np.sum(np.multiply(a[i:i + vh, j:j + vw], v))
          for j in range(aw - vw + 1)] for i in range(ah - vh + 1)]
    return k

def aconv(a, v):
    n, m, h, w = np.shape(a)
    vn, vm, vh, vw = np.shape(v)
    d, s = h -vh +1, w -vw +1
    z = np.zeros((n, vn, d, s))
    for i in range(n):
        for j in range(vn):
            for k in range(m):
                    z[i][j] += conv2d(a[i][k], v[j][k])
    return z


def sconv(a, v):
    n, m, h, w = np.shape(a)
    vn, vm, vh, vw = np.shape(v)
    d, s = h -vh +1, w -vw +1
    z = np.zeros((n, vn, d, s))
    for i in range(n):
        for j in range(vn):
            for k in range(m):
                    z[i][j] += sg.convolve2d(a[i][k], np.rot90(v[j][k],2), mode='valid')
    return z

def agconv():
    return aconv(a,b)

def sgconv():
    return sconv(a, b)

def imconv():
    return imcol(a, b)

if __name__=='__main__':
    print sgconv()[0, 0, 0]
    print imconv()[0, 0, 0]
    print timeit.timeit('agconv()',setup='from __main__ import agconv',number=10)
    print timeit.timeit('sgconv()',setup='from __main__ import sgconv',number=10)
    print timeit.timeit('imconv()', setup='from __main__ import imconv', number=10)

"""
[ 8.5865   8.59145  8.5964   8.60135  8.6063   8.61125  8.6162   8.62115
  8.6261   8.63105  8.636    8.64095  8.6459   8.65085  8.6558   8.66075
  8.6657   8.67065  8.6756   8.68055  8.6855   8.69045  8.6954   8.70035]
[ 8.5865   8.59145  8.5964   8.60135  8.6063   8.61125  8.6162   8.62115
  8.6261   8.63105  8.636    8.64095  8.6459   8.65085  8.6558   8.66075
  8.6657   8.67065  8.6756   8.68055  8.6855   8.69045  8.6954   8.70035]
24.7670379661
0.686509787089
0.523622603254

"""

在反向传播的时候需要一个相应的col2im。

import numpy as np

a=range(3*2*20*20)
a=np.reshape(a,(3,2,20,20))*0.001
image_shape=a.shape
filter_shape=[2,2,5,5]
def im2col(a):
    n, m, h, w = image_shape
    vn, vm, vh, vw = filter_shape
    t = h - vh + 1
    kh = t ** 2
    kw = vh * vh
    z = np.zeros((n, kh, kw * m))
    for i in range(n):
        for j in range(m):
            kt = []
            for it in range(t):
                for jt in range(t):
                    kt.append((a[i, j, it:it + vh, jt:jt + vh]).flatten())
            z[i, :, j * kw:j * kw + kw] = kt
    z = np.reshape(z, (n * kh, kw * m))
    return z

def col2im(a):
    n, m, h, w = image_shape
    vn, vm, vh, vw = filter_shape
    t = h - vh + 1
    kh = t ** 2
    kw = vh * vh
    a = np.reshape(a, (n, kh, kw * m))
    z=np.zeros((n,m,h,h))
    for i in range(n):
        for k in range(m):
            for j in range(h):
                rh = vh * min(j, vh - 1)
                b = max(j + 1 - vh, 0) * t
                z[i,k,j] = np.append(a[i,b,rh:rh + vh], a[i, b + 1:b + t, rh + vh - 1])
    return z

x=im2col(a)
t=col2im(x)
print t[0][0]==a[0,0]

"""
[[ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]
 [ True  True  True  True  True  True  True  True  True  True  True  True
   True  True  True  True  True  True  True  True]]
"""

完整的代码如下。

# coding:utf8
import cPickle
import numpy as np

class ConvPoolLayer(object):
    def __init__(self, image_shape,filter_shape,poolsize=(2,2)):
        self.filter_shape = filter_shape
        self.image_shape = image_shape
        self.w = np.random.normal(loc=0, scale=np.sqrt(1.0/np.prod(filter_shape[1:])),
                                  size=(np.prod(filter_shape[1:]),filter_shape[0]))
        self.b = np.random.normal(loc=0, scale=0.2, size=(1,filter_shape[0],1,1))
        self.poolsize = poolsize
        self.samp_shape=(image_shape[-2] - filter_shape[-2] + 1,image_shape[-1] - filter_shape[-1] + 1)
        self.out_shape=(self.samp_shape[0]/poolsize[0],self.samp_shape[1]/poolsize[1])

    def im2col(self,a):
        n, m, h, w = self.image_shape
        vn, vm, vh, vw = self.filter_shape
        t = h - vh + 1
        kh = t ** 2
        kw = vh * vh
        z = np.zeros((n, kh, kw * m))
        for i in range(n):
            for j in range(m):
                kt = []
                for it in range(t):
                    for jt in range(t):
                        kt.append((a[i, j, it:it + vh, jt:jt + vh]).flatten())
                z[i, :, j * kw:j * kw + kw] = kt
        z = np.reshape(z, (n * kh, kw * m))
        return z

    def col2im(self,a):
        n, m, h, w = self.image_shape
        vn, vm, vh, vw = self.filter_shape
        t = h - vh + 1
        kh = t ** 2
        kw = vh * vh
        a = np.reshape(a, (n, kh, kw * m))
        z = np.zeros((n, m, h, h))
        for i in range(n):
            for k in range(m):
                for j in range(h):
                    rh = vh * min(j, vh - 1)
                    b = max(j + 1 - vh, 0) * t
                    z[i, k, j] = np.append(a[i, b, rh:rh + vh], a[i, b + 1:b + t, rh + vh - 1])
        return z

    def fw_shape(self,a):
        res = np.reshape(a, (self.image_shape[0], -1, self.filter_shape[0]))
        res = np.rollaxis(res, 2, 1)
        res = np.reshape(res, (self.image_shape[0], -1, self.samp_shape[0], self.samp_shape[1]))
        return res

    def bp_shape(self,a):
        res = np.reshape(a, (self.image_shape[0], self.filter_shape[0], -1))
        res = np.rollaxis(res, 1, 3)
        res = np.reshape(res, (-1, self.filter_shape[0]))
        return res

    def feedforward(self, a):
        z = self.im2col(a)
        res = np.dot(z, self.w)
        res = self.fw_shape(res)
        self.out = self.relu(res+self.b)
        return np.array(self.pool2d(self.out))

    def backprop(self, x, dnext,eta=0.001):
        if dnext.ndim<3:
            dnext = np.reshape(dnext,(self.image_shape[0],self.filter_shape[0], self.out_shape[0], self.out_shape[1]))
        u = self.relu_prime(self.out)
        dnext = np.multiply(u,self.up(dnext,self.poolsize[0]))
        delta = self.bp_shape(dnext)/self.image_shape[0]
        x=self.im2col(x)
        out_delta = np.dot(delta, self.w.T)
        out_delta=self.col2im(out_delta)
        w = np.dot(x.T, delta)
        self.w -= eta * w
        self.b -= eta * np.reshape(np.sum(delta,0),(self.b.shape))
        return out_delta

    def pool2d(self,input, ds=(2, 2), mode='max'):
        fun = np.max
        if mode == 'sum':
            fun = np.sum
        elif mode == 'average':
            fun = np.average
        n, m, h, w = np.shape(input)
        d, s = ds
        zh = h / d + h % d
        zw = w / s + w % s
        z = np.zeros((n, m, zh, zw))
        for k in range(n):
            for o in range(m):
                for i in range(zh):
                    for j in range(zw):
                        z[k, o, i, j] = fun(input[k, o, d * i:min(d * i + d, h), s * j:min(s * j + s, w)])
        return z

    def up(self,a,n):
        b=np.ones((n,n))
        return np.kron(a,b)

    def relu(self,z):
        return np.maximum(z, 0.0)

    def relu_prime(self,z):
        z[z>0]=1
        return z

class SoftmaxLayer(object):
    def __init__(self, in_num=100,out_num=10):
        self.weights = np.random.randn(in_num, out_num)/np.sqrt(out_num)

    def feedforward(self, input):
        self.out=self.softmax(np.dot(input, self.weights))
        return self.out

    def backprop(self, input, y,eta=0.001):
        o=self.out
        delta =o-y
        out_delta=np.dot(delta,self.weights.T)
        w = np.dot(input.T,delta)
        self.weights-=eta*(w)
        return out_delta

    def softmax(self,a):
        m = np.exp(a)
        return m / (np.sum(m,axis=1)[:,np.newaxis])

class FullLayer(object):
    def __init__(self, in_num=720,out_num=100):
        self.in_num=in_num
        self.out_num=out_num
        self.biases = np.random.randn(out_num)
        self.weights = np.random.randn(in_num, out_num)/np.sqrt(out_num)

    def feedforward(self, x):
        if x.ndim>2:
            x = np.reshape(x, (len(x), self.in_num))
        self.out = self.sigmoid(np.dot(x, self.weights)+self.biases)
        return self.out

    def backprop(self, x,delta,eta=0.001):
        if x.ndim>2:
            x = np.reshape(x, (len(x), self.in_num))
        sp=self.sigmoid_prime(self.out)
        delta = delta * sp
        out_delta = np.dot(delta, self.weights.T)
        w = np.dot( x.T,delta)
        self.weights-=eta*w
        self.biases -= eta*np.sum(delta,0)
        return out_delta

    def sigmoid(self,z):
        return 1.0/(1.0+np.exp(-z))

    def sigmoid_prime(self,z):
        return z*(1-z)

class Network(object):
    def __init__(self, layers):
        self.layers=layers
        self.num_layers = len(layers)
        self.a=[]

    def feedforward(self, x):
        self.a.append(x)
        for layer in self.layers:
            x=layer.feedforward(x)
            self.a.append(x)
        return x

    def SGD(self, training_data, test_data,epochs, mini_batch_size, eta=0.001):
        self.n = len(training_data[0])
        self.mini_batch_size=mini_batch_size
        self.eta = eta
        cx=range(epochs)
        for j in cx:
            for k in xrange(0, self.n , mini_batch_size):
                batch_x = np.array(training_data[0][k:k + mini_batch_size])
                batch_y = training_data[1][k:k + mini_batch_size]
                self.backprop(batch_x,batch_y)
                if k%500==0:
                    print "Epoch {0}:{1} cost={3}, test:{2}".format(j,k,
                    self.evaluate([test_data[0],test_data[1]]),self.cost)

    def backprop(self, x_in, y):
        self.feedforward(x_in)
        for i in range(self.num_layers):
            delta=self.layers[-i-1].backprop(self.a[-i-2],y,eta=self.eta)
            y=delta

    def evaluate(self, test_data):
        x,y=test_data
        num=len(x)
        size=self.mini_batch_size
        x=[self.feedforward(np.array(x[size*i:size*i+size])) for i in range((num/size))]
        x=np.reshape(x,(num,np.shape(x)[-1]))
        xp = np.argmax(x, axis=1)
        yp= np.argmax(y, axis=1) if y[0].ndim else y
        self.cost = -np.mean(np.log(x)[np.arange(num),yp])
        return np.mean(yp == xp)*100


if __name__ == '__main__':
        def get_data(data):
            return [np.reshape(x, (1,28,28)) for x in data[0]]

        def get_label(i):
            c = np.zeros((10))
            c[i] = 1
            return c

        f = open('data/mnist.pkl', 'rb')
        training_data, validation_data, test_data = cPickle.load(f)
        training_inputs = get_data(training_data)
        training_label=[get_label(y_) for y_ in training_data[1]]
        test_inputs = get_data(test_data)
        test = zip(test_inputs,test_data[1])
        size=10
        net = Network([ConvPoolLayer(image_shape=[size,1,28,28],filter_shape=[8,1,5,5],poolsize=(2,2)),
                       ConvPoolLayer(image_shape=[size,8,12,12],filter_shape=[16,8,5,5], poolsize=(2,2)),
                       FullLayer(in_num=16*4*4,out_num=100),
                       SoftmaxLayer(in_num=100,out_num=10)])
        net.SGD([training_inputs,training_label],[test_inputs[:500],test_data[1][:500]],
                epochs=3,mini_batch_size=size, eta=0.05)

        # Epoch 0:14000 cost=0.134660777083, test:96.0