【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得到的全部結果相加后得到最后的

for i = 1:numImages
    for fil2 = 1:numFilters2
        convolvedImage = zeros(convDim, convDim);
        for fil1 = 1:numFilters1
            filter = squeeze(W(:,:,fil1,fil2));
            filter = rot90(squeeze(filter),2);
            im = squeeze(images(:,:,fil1,i));
            convolvedImage = convolvedImage + conv2(im,filter,'valid');
       end
        convolvedImage = bsxfun(@plus,convolvedImage,b(fil2));
        convolvedImage = 1 ./ (1+exp(-convolvedImage));
        convolvedFeatures(:, :, fil2, i) = convolvedImage;
     end
end      

最后整個Feedforward的過程代碼如下：

%Feedfoward Pass
activations1 = cnnConvolve4D(mb_images, Wc1, bc1);
activationsPooled1 = cnnPool(poolDim1, activations1);
activations2 = cnnConvolve4D(activationsPooled1, Wc2, bc2);
activationsPooled2 = cnnPool(poolDim2, activations2);

% Reshape activations into 2-d matrix, hiddenSize x numImages,
% for Softmax layer
activationsPooled2 = reshape(activationsPooled2,[],numImages);

%% --------- Softmax Layer ---------
probs = exp(bsxfun(@plus, Wd * activationsPooled2, bd));
sumProbs = sum(probs, 1);
probs = bsxfun(@times, probs, 1 ./ sumProbs);
            

3. Caculate cost

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

logp = log(probs);
index = sub2ind(size(logp),mb_labels',1:size(probs,2));
ceCost = -sum(logp(index));        
wCost = lambda/2 * (sum(Wd(:).^2)+sum(Wc1(:).^2)+sum(Wc2(:).^2));
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和我們所預測的結果之間的差值：

output = zeros(size(probs));
output(index) = 1;
DeltaSoftmax = (probs - output);
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函數共同實現，具體的解釋參考這里。代碼如下：

DeltaUnpool2 = zeros(convDim2,convDim2,numFilters2,numImages);        
for imNum = 1:numImages
    for FilterNum = 1:numFilters2
        unpool = DeltaPool2(:,:,FilterNum,imNum);
        DeltaUnpool2(:,:,FilterNum,imNum) = kron(unpool,ones(poolDim2))./(poolDim2 ^ 2);
    end
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了。如下圖所示：

代碼如下：

%error of first pooling layer
DeltaPooled1 = zeros(outputDim1,outputDim1,numFilters1,numImages);
for i = 1:numImages
   for f1 = 1:numFilters1
      for f2 = 1:numFilters2
         DeltaPooled1(:,:,f1,i) = DeltaPooled1(:,:,f1,i) + convn(DeltaConv2(:,:,f2,i),Wc2(:,:,f1,f2),'full');
      end
   end
end

然后同樣進行上面解釋過的unpool過程得到DeltaUnpooled1:

%error of first convolutional layer
DeltaUnpool1 = zeros(convDim1,convDim1,numFilters1,numImages);
for imNum = 1:numImages
     for FilterNum = 1:numFilters1
         unpool = DeltaPooled1(:,:,FilterNum,imNum);
         DeltaUnpool1(:,:,FilterNum,imNum) = kron(unpool,ones(poolDim1))./(poolDim1 ^ 2);
      end
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計算就可以了，代碼如下：

 % softmax layer
 Wd_grad = DeltaSoftmax*activationsPooled2';
 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對應的梯度了，如下圖所示：

代碼如下：

for fil2 = 1:numFilters2
    for fil1 = 1:numFilters1
        for im = 1:numImages
            Wc2_grad(:,:,fil1,fil2) = Wc2_grad(:,:,fil1,fil2) + conv2(activationsPooled1(:,:,fil1,im),rot90(DeltaConv2(:,:,fil2,im),2),'valid');
        end
    end
    temp = DeltaConv2(:,:,fil2,:);
    bc2_grad(fil2) = sum(temp(:));
end

4.2.3 Wc1，bc1的梯度計算。

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

 % first convolutional layer
 for fil1 = 1:numFilters1
     for channel = 1:imageChannel
         for im = 1:numImages
             Wc1_grad(:,:,channel,fil1) = Wc1_grad(:,:,channel,fil1) + conv2(mb_images(:,:,channel,im),rot90(DeltaConv1(:,:,fil1,im),2),'valid');
         end
     end
     temp = DeltaConv1(:,:,fil1,:);
     bc1_grad(fil1) = sum(temp(:));
 end

4.2.4 梯度值更新。

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

 Wd_velocity = mom*Wd_velocity + alpha*(Wd_grad/minibatch+lambda*Wd);
 bd_velocity = mom*bd_velocity + alpha*bd_grad/minibatch;
 Wc2_velocity = mom*Wc2_velocity + alpha*(Wc2_grad/minibatch+lambda*Wc2);
 bc2_velocity = mom*bc2_velocity + alpha*bc2_grad/minibatch;
 Wc1_velocity = mom*Wc1_velocity + alpha*(Wc1_grad/minibatch+lambda*Wc1);
 bc1_velocity = mom*bc1_velocity + alpha*bc1_grad/minibatch;
                        
 Wd = Wd - Wd_velocity;
 bd = bd - bd_velocity;
 Wc2 = Wc2 - Wc2_velocity;
 bc2 = bc2 - bc2_velocity;
 Wc1 = Wc1 - Wc1_velocity;
 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