【Deep Learning】兩層CNN的MATLAB實現

想本身動手寫一個CNN好久了，論文和代碼之間的差距有一個銀河系那麼大。

在實現兩層的CNN以前，首先實現了UFLDL中與CNN有關的做業。而後參考它的代碼搭建了一個一層的CNN。最後實現了一個兩層的CNN，碼代碼花了一天，調試花了5天，我也是醉了。這裏記錄一下經過代碼對CNN加深的理解。

首先，dataset是MNIST。這裏層的概念是指convolution+pooling，有些地方會把convolution和pooling分別做爲兩層看待。

1.CNN的結構

這個兩層CNN的結構以下：

圖一

各個變量的含義以下（和代碼中的變量名是一致的）

images:輸入的圖片，一張圖片是28*28，minibatch的大小設置的是150，因此輸入就是一個28*28*150的矩陣。

Wc1,bc1：第一層卷積的權重和偏置。一共8個filter，每一個大小爲5*5。

activations1：經過第一層卷積獲得的feature map，大小爲（28-5+1）*（28-5+1）*8*150，其中8是第一層卷積filter的個數，150是輸入的image的個數。

activationsPooled1：將卷積後的feature map進行採樣後的feature map，大小爲(24/2)*(24/2)*8*150。

Wc2,bc2：第二層卷積的權重和偏置。一共10個filter，每一個大小爲5*5*8.注意第二層的權重是三維的，這就是兩層卷積網絡和一層卷積網絡的差異，對於每張圖像，第一層輸出的對應這張圖像的feature map是12*12*8,而第二層的一個filter和這個feature map卷積後獲得一張8*8的feature map，因此第二層的filter都是三維的。具體操做後面再詳細介紹。

activations2：經過第二層卷積獲得的feature map，大小爲（12-5+1）*（12-5+1）*10*150，其中10是第二層卷積filter的個數，150是輸入的image的個數。

activationsPooled2：將卷積後的feature map進行採樣後的feature map，大小爲(8/2)*(8/2)*10*150。

activationsPooled2'：第二層卷積完了以後，要把一張image對應的全部feature map reshape成一列，那麼這一列的長度就是4*4*10=160，因此reshape後獲得一個160*150的大feature map.（代碼中仍然是用的activationsPooled2）。

Wd,bd：softmax層的權重和偏執。

probs：對全部圖像所屬分類的預測結果，每一列對應一張圖像，一共10行，第i行表明這張圖像屬於第i類的機率。

---

 

從實現的角度來講，一個CNN主要能夠分紅三大塊：Feedfoward Pass，Caculate cost和 Backpropagation.這裏就詳細介紹這三塊。

2. Feedfoward Pass

這個過程主要是輸入一張圖像，經過目前的權重，使得圖像依次經過每一層的convolution或者pooling操做，最後獲得對圖像分類的機率預測。這裏比較tricky的部分是對三維的feature map的卷積過程。好比說對於第一層pooling輸出的feature map，一張圖像對應一個尺寸爲12*12*8的feature map，這個時候要用第二層卷積層中一個5*5*8的filter(Wc2中的一個filter)和它進行卷積，最終獲得一個8*8的feature map。整個過程能夠用圖二，圖三來表示：

　　　　　　　　　　　　　　　　圖二

繼續放大，來看看輸入的feature map是怎樣和第二層的一個filter進行卷積的：

　　　　　　　　　　　　　　圖三

也就是說，這個卷積是經過filter中每個filter:Wc2(:,:,fil1,fil2)和activationsPooled1(:,:,fil1,imageNum)卷積，而後將全部的fil1=1:8獲得的所有結果相加後獲得最後的

 1 for i = 1:numImages
 2     for fil2 = 1:numFilters2
 3         convolvedImage = zeros(convDim, convDim);
 4         for fil1 = 1:numFilters1
 5             filter = squeeze(W(:,:,fil1,fil2));
 6             filter = rot90(squeeze(filter),2);
 7             im = squeeze(images(:,:,fil1,i));
 8             convolvedImage = convolvedImage + conv2(im,filter,'valid');
 9        end
10         convolvedImage = bsxfun(@plus,convolvedImage,b(fil2));
11         convolvedImage = 1 ./ (1+exp(-convolvedImage));
12         convolvedFeatures(:, :, fil2, i) = convolvedImage;
13      end
14 end      

最後整個Feedforward的過程代碼以下：

 1 %Feedfoward Pass
 2 activations1 = cnnConvolve4D(mb_images, Wc1, bc1);
 3 activationsPooled1 = cnnPool(poolDim1, activations1);
 4 activations2 = cnnConvolve4D(activationsPooled1, Wc2, bc2);
 5 activationsPooled2 = cnnPool(poolDim2, activations2);
 6 
 7 % Reshape activations into 2-d matrix, hiddenSize x numImages,
 8 % for Softmax layer
 9 activationsPooled2 = reshape(activationsPooled2,[],numImages);
10 
11 %% --------- Softmax Layer ---------
12 probs = exp(bsxfun(@plus, Wd * activationsPooled2, bd));
13 sumProbs = sum(probs, 1);
14 probs = bsxfun(@times, probs, 1 ./ sumProbs);
15             

3. Caculate cost

這裏使用的loss function是cross entropy function，關於這個函數的細節能夠看這裏。這個函數比起squared error function的好處是它在表現越差的時候學習的越快。這和咱們的直覺是相符的。而對於squared error function，它在loss比較大的時候反而進行的梯度更新值很小，即學習的很慢，具體解釋也參見上述連接。計算cost的代碼十分直接，這裏直接使用了ufldl做業中的代碼。

1 logp = log(probs);
2 index = sub2ind(size(logp),mb_labels',1:size(probs,2));
3 ceCost = -sum(logp(index));        
4 wCost = lambda/2 * (sum(Wd(:).^2)+sum(Wc1(:).^2)+sum(Wc2(:).^2));
5 cost = ceCost/numImages + wCost;

注意ceCost是loss function真正的cost，而wCost是weight decay引發的cost，咱們指望學習到的網絡的權重都偏小，對於這一點如今沒有很完備的解釋，咱們指望權值比較小的一個緣由是小的權值使得輸入波動比較大的時候，網絡的各部分的值變化不至於太大，不然網絡會不穩定。

4. Backpropagation

Backpropagation算法其實能夠分紅兩部分：計算error和gradient

4.1 caculate error of each layer

對於每一層生成的feature map，咱們計算一個偏差（殘差），說明這一層計算出來的結果和它應該給出的「正確」結果之間的差值。

計算這個偏差的過程，其實就是「找源頭」的過程，只要知道某張feature map和該層的哪些filter生成了下一層的哪些feature map，而後用下一層feature map對應的偏差和filter就能夠獲得要計算的feature map對應的偏差了。

那麼對於最後一層softmax的偏差就很好理解了，就是ground truth的labels和咱們所預測的結果之間的差值：

1 output = zeros(size(probs));
2 output(index) = 1;
3 DeltaSoftmax = (probs - output);
4 t = -DeltaSoftmax;

output是把ground truth的labels整成一個10*150的矩陣，output(i,j)=1表示圖像j屬於第i類，output(i,j)=0表示圖像j不屬於第i類。

接下來把這個殘差一層層的依次推回到pooling2->convolution2->pooling1->convolution1這些層。

用到的公式是如下四個，具體的推倒參見這裏，下面也會有一部分對這些公式直觀的解釋。

4.1.1 pooling2層的偏差

這一層比較簡單，根據上面的公式BP2，咱們直接能夠用Wd' * DeltaSoftmax就能夠獲得這一層的偏差(由於這一層沒有sigmoid函數，因此沒有BP2後面的導數部分)。固然，Wd是一個10*160的矩陣，而DeltaSoftmax和probs同維度，即10*150（參見圖一），它們相乘後獲得160*150的矩陣，其中每一列對應一張圖像的偏差。而咱們所要求得pooling2層的feature map的維度是4*4*10*150，因此咱們要把獲得的160*150的矩陣reshape成pooling2的feature map所對應的偏差。具體代碼以下：

DeltaPool2 = reshape(Wd' * DeltaSoftmax,outputDim2,outputDim2,numFilters2,numImages);

獲得pooling2層的偏差後，一個很重要的操做是unpool的過程。爲何要用這個過程呢？咱們先來看一個簡單的pooling過程：

假設有一張4*4的feature map，對它進行average pooling：

在上面的pooling過程當中，採樣後的featuremap中紅色部分的值來自於未採樣的feature map中紅色部分的4個值取平均，因此紅色部分的值的偏差，就由這4個紅色的值「負責」，這個偏差在unpool的過程當中就在這4個值對應的error上均分。其餘顏色的部分一樣的道理。unpool部分的代碼由conv函數和kron函數共同實現，具體的解釋參考這裏。代碼以下：

1 DeltaUnpool2 = zeros(convDim2,convDim2,numFilters2,numImages);        
2 for imNum = 1:numImages
3     for FilterNum = 1:numFilters2
4         unpool = DeltaPool2(:,:,FilterNum,imNum);
5         DeltaUnpool2(:,:,FilterNum,imNum) = kron(unpool,ones(poolDim2))./(poolDim2 ^ 2);
6     end
7 end

4.1.2 convolution2的偏差

有了上述unpool的偏差，咱們就能夠直接用公式BP2計算了：

DeltaConv2 = DeltaUnpool2 .* activations2 .* (1 - activations2);

其中的activations2 .* (1 - activations2)對應BP2中的σ'(z)，這是sigmoid函數一個很好的性質。

4.1.3 pooling1的偏差

這一層的偏差計算的關鍵一樣是找準「源頭」。在feedfoward的過程當中，第二層的convolution過程以下：

假設咱們要求第一張圖像的第二張feature map（黑色那張）對應的偏差error。那麼咱們只要搞清楚它「幹了多少壞事」，而後把這些「壞事」加起來就是它的偏差了。如上圖所示，假設第二層convolution由4個5*5*3的filter，那麼這張黑色的feature map分別和這4個filter中每一個filters的第二張filter進行過卷積，而且這些卷積的結果分別貢獻給了第一張圖convolved feature map的第1，2，3，4個feature map（上圖convolved feature map中和filter顏色對應的那幾張feature map）。因此，要求紫色的feature map的偏差，咱們用convolved feature map中對應顏色的error和對應顏色filter卷積，而後將全部的卷積結果相加，就能夠獲得紫色的feature map的error了。以下圖所示：

代碼以下：

1 %error of first pooling layer
2 DeltaPooled1 = zeros(outputDim1,outputDim1,numFilters1,numImages);
3 for i = 1:numImages
4    for f1 = 1:numFilters1
5       for f2 = 1:numFilters2
6          DeltaPooled1(:,:,f1,i) = DeltaPooled1(:,:,f1,i) + convn(DeltaConv2(:,:,f2,i),Wc2(:,:,f1,f2),'full');
7       end
8    end
9 end

而後一樣進行上面解釋過的unpool過程獲得DeltaUnpooled1:

1 %error of first convolutional layer
2 DeltaUnpool1 = zeros(convDim1,convDim1,numFilters1,numImages);
3 for imNum = 1:numImages
4      for FilterNum = 1:numFilters1
5          unpool = DeltaPooled1(:,:,FilterNum,imNum);
6          DeltaUnpool1(:,:,FilterNum,imNum) = kron(unpool,ones(poolDim1))./(poolDim1 ^ 2);
7       end
8 end

4.1.4 convolution1的偏差

這層的偏差仍是根據BP2公式計算：

DeltaConv1 = DeltaUnpool1 .* activations1 .* (1-activations1);

4.2 Caculate Gradient

對於梯度的計算相似於偏差的計算，對於每個filter，找到它爲哪些feature map的計算「作出貢獻」，而後用這些feature map的偏差計算相應的梯度並求和。在咱們的CNN中，有三個權值Wc1,Wc2,Wd,bc1,bc2,bd的梯度須要計算。

4.2.1 Wd，bd的梯度計算。

這兩個梯度十分好計算，只要根據公式BP4計算就能夠了，代碼以下：

1  % softmax layer
2  Wd_grad = DeltaSoftmax*activationsPooled2';
3  bd_grad = sum(DeltaSoftmax,2);

4.2.2 Wc2，bc2的梯度計算。

計算Wc2，首先要知道在forward pass過程當中，Wc2中的某個filter生成了哪些feature map，而後在用這些feature map的error來計算filter的梯度，feedforward的過程以下圖所示：

假設咱們要計算黑色的filter對應的梯度，在feedforward的過程當中，這個filter和左邊的pooling1層輸出的feature map卷積，生成右邊對應顏色的feature map，那麼在backpropagation的過程當中，咱們就用右邊這些feature map對應的偏差error和左邊輸入的feature map卷積，最後把每張圖像的卷積結果相加，就能夠獲得黑色的filter對應的梯度了，以下圖所示：

代碼以下：

1 for fil2 = 1:numFilters2
2     for fil1 = 1:numFilters1
3         for im = 1:numImages
4             Wc2_grad(:,:,fil1,fil2) = Wc2_grad(:,:,fil1,fil2) + conv2(activationsPooled1(:,:,fil1,im),rot90(DeltaConv2(:,:,fil2,im),2),'valid');
5         end
6     end
7     temp = DeltaConv2(:,:,fil2,:);
8     bc2_grad(fil2) = sum(temp(:));
9 end

4.2.3 Wc1，bc1的梯度計算。

Wc1和bc1的計算和Wc2，bc2是相似的，只是第一層的輸入是圖像數據集images，因此只要把上述代碼中的numFilters2換成numFilters1,numFilters1換成imageChannel(圖像的通道RGB圖像對應3，這裏用灰度圖像，因此imageChannel的值是1)，activationsPooled1換成圖像集mb_images，DeltaConv2換成DeltaConv1就能夠了。

 1  % first convolutional layer
 2  for fil1 = 1:numFilters1
 3      for channel = 1:imageChannel
 4          for im = 1:numImages
 5              Wc1_grad(:,:,channel,fil1) = Wc1_grad(:,:,channel,fil1) + conv2(mb_images(:,:,channel,im),rot90(DeltaConv1(:,:,fil1,im),2),'valid');
 6          end
 7      end
 8      temp = DeltaConv1(:,:,fil1,:);
 9      bc1_grad(fil1) = sum(temp(:));
10  end

4.2.4 梯度值更新。

這一步就十分簡單了，直接用計算好的梯度去更新權重就能夠了。不過使用了衝量和weight decay，其中alpha是學習速率，lambda是weight decay factor，代碼以下：

 1  Wd_velocity = mom*Wd_velocity + alpha*(Wd_grad/minibatch+lambda*Wd);
 2  bd_velocity = mom*bd_velocity + alpha*bd_grad/minibatch;
 3  Wc2_velocity = mom*Wc2_velocity + alpha*(Wc2_grad/minibatch+lambda*Wc2);
 4  bc2_velocity = mom*bc2_velocity + alpha*bc2_grad/minibatch;
 5  Wc1_velocity = mom*Wc1_velocity + alpha*(Wc1_grad/minibatch+lambda*Wc1);
 6  bc1_velocity = mom*bc1_velocity + alpha*bc1_grad/minibatch;
 7                         
 8  Wd = Wd - Wd_velocity;
 9  bd = bd - bd_velocity;
10  Wc2 = Wc2 - Wc2_velocity;
11  bc2 = bc2 - bc2_velocity;
12  Wc1 = Wc1 - Wc1_velocity;
13  bc1 = bc1 - bc1_velocity;

總結如下Backpropagation算法的關鍵，就是找準偏差的源頭，而後將這個偏差順着源頭從最後一層推回到第一層，沿路根據偏差更新權重，以此來訓練神經網絡。

5. CNN效果及代碼

如下是cost在迭代過程當中的變化圖像，一共進行了3次迭代，每次對400個minibatch進行stochastic gradient descent，每一個minibatch有150張圖像。

最後的結果：

代碼參見個人github。

有趣的一點是，以前我用一層的CNN在MNIST上能夠達到97.4%的準確率，換成這個兩層的CNN，準確率卻降低了。一種多是我沒有進行fine tuning，上述的參數有些是參考別人的，有些是本身隨便設置的；另外一個緣由多是overfitting，在沒有加入weight decay的代碼以前，獲得的準確率只有94%，weight decay減輕了部分overfitting的現象，準確率提升了兩個百分點。

以上是我的實現CNN的筆記，歡迎大神指正。

參考資料：

【1】http://neuralnetworksanddeeplearning.com/

【2】https://github.com/yaolubrain/cnn_linear_max

【3】https://github.com/jakezhaojb/toy-cnn/tree/master/matlab