這個exercise需要完成cnn中的forward pass,cost,error和gradient的計算。需要弄清楚每一層的以上四個步驟的原理，並且要充分利用matlab的矩陣運算。大概把過程總結了一下如下圖所示：

STEP 1：Implement CNN Objective

STEP 1a: Forward Propagation

Forward Propagation主要是為了計算輸入圖片經過神經網路後的輸出，這個網路有三層：convolution->pooling->softmax(dense connected)，卷積層對於每個影象用所有的模板進行卷積；pooling層對卷積層輸出的feature map取樣；softmax層根據pooling層輸出的feature預測影象的分類結果。其中convolution和pooling操作在

%%% YOUR CODE HERE %%%
%呼叫之前已經實現的函式
activations = cnnConvolve(filterDim, numFilters, images, Wc, bc);%sigmoid(wx+b)
activationsPooled = cnnPool(poolDim, activations);

% Reshape activations into 2-d matrix, hiddenSize x numImages,
% for Softmax layer
 
%將activationsPooled從outDim*outDim*numFilters*numImages拼接成hiddenSize*numImages的大矩陣
activationsPooled = reshape(activationsPooled,[],numImages);

%% Softmax Layer
%  Forward propagate the pooled activations calculated above into a
%  standard softmax layer. For your convenience we have reshaped
 
%  activationPooled into a hiddenSize x numImages matrix.  Store the
%  results in probs.

% numClasses x numImages for storing probability that each image belongs to
% each class.
probs = zeros(numClasses,numImages);

%%% YOUR CODE HERE %%%
h = exp(bsxfun(@plus,Wd * activationsPooled,bd));
probs = bsxfun(@rdivide,h,sum(h,1));

STEP 1b: Calculate Cost

計算梯度下降要優化的目標函式，主要分為兩部分，一部分是由於分類器輸出結果和真實結果的差異引起的誤差函式，另一部分是對權重w的正則約束。第一部分可以參考softmax regression中對損失函式的計算，第二部分就是對Wc和Wd的所有項求平方和。類似下面的公式，不過第一項中的J是softmax的cross entropy損失函式。最後要對第一項除以影象的總數，這是十分重要的，一開始我沒有除，最後得到的演算法是發散的，原因可能是第一項數值過大，直接把正則項的影響給忽略了。

程式碼：

%%% YOUR CODE HERE %%%
logp = log(probs);
index = sub2ind(size(logp),labels',1:size(probs,2));
ceCost = -sum(logp(index));
wCost = lambda/2 * (sum(Wd(:).^2)+sum(Wc(:).^2));
cost = ceCost/numImages + wCost;

STEP 1c: Backpropagation

BP演算法首先要計算各層的對最終誤差的貢獻delta。

softmax層：這一層的誤差最容易計算，只要用ground truth減去神經網路的輸出probs就可以了：

output = zeros(size(probs));
output(index) = 1;
DeltaSoftmax = probs - output;

pool層：這一層首先根據公式δ^l ＝ Wδ^l+1 * f'(z^l)（pool層沒有f'(z^l)這一項）計算該層的error，此時得到一個hiddenSize*numImages的矩陣，首先利用reshape函式把error還原成一個convDim*convDim*numFilters*numImages的矩陣，在pooling操作時，pooling層一個節點的輸入是conv層2*2個節點的輸出（假設poolDim=2）如下圖所示：

所以pooling層的這個節點要將自己的error在這2*2個節點中平均分（因為使用的是mean pooling）,UFLDL上面提示了可以用kron這個函式來實現，這樣如上圖所示，就可以通過pooling層一個2*2的filter對應的error計算得到convolution層一個4*4的filter對應的error了。程式碼如下：

DeltaPool = reshape(Wd' * DeltaSoftmax,outputDim,outputDim,numFilters,numImages);
DeltaUnpool = zeros(convDim,convDim,numFilters,numImages);

for imNum = 1:numImages
    for FilterNum = 1:numFilters
        unpool = DeltaPool(:,:,FilterNum,imNum);
        DeltaUnpool(:,:,FilterNum,imNum) = kron(unpool,ones(poolDim))./(poolDim ^ 2);
    end
end

convolution層：還是根據公式δ^l ＝ Wδ^l+1 * f'(z^l)來計算：

DeltaConv = DeltaUnpool .* activations .* (1 - activations);

STEP 1d: Gradient Calculation

整個cnn一共有三層：convolution->pooling->softmax(dense connected)，只有convolution和softmax層有權重，分別為Wc，bc，Wd，bd。那麼就要計算目標函式J對他們的倒數以便在梯度下降中更新W和b。

Wd和bd的梯度計算：

根據下面兩個公式：

其中a^l-1對應pooling層的激勵（輸出）activitonsPooled，δ^l就是這一層的誤差DeltaSoftmax，程式碼如下：

Wd_grad = (1./numImages) .* DeltaSoftmax*activationsPooled'+lambda*Wd;
bd_grad = (1./numImages) .* sum(DeltaSoftmax,2);

Wc和bc的梯度計算：

還是根據上面兩個計算梯度的公式，不過麻煩就麻煩在l-1層其實是輸入的影象，所以a^l-1對應的是輸入的影象，那麼就得用for迴圈逐個便利影象並利用UFLDL上提供的公式計算對應梯度：

這裡為了方便，先對所有DeltaConv進行旋轉，然後再用for迴圈依次求出梯度：

%%% YOUR CODE HERE %%%
Wd_grad = (1./numImages) .* DeltaSoftmax*activationsPooled'+lambda*Wd;
bd_grad = (1./numImages) .* sum(DeltaSoftmax,2);

bc_grad = zeros(size(bc));
Wc_grad = zeros(filterDim,filterDim,numFilters);

for filterNum = 1:numFilters
    error = DeltaConv(:,:,filterNum,:);
    bc_grad(filterNum) = (1./numImages) .* sum(error(:));
end

%旋轉所有DealtaConv
for filterNum = 1:numFilters
    for imNum = 1:numImages
        error = DeltaConv(:,:,filterNum,imNum);
        DeltaConv(:,:,filterNum,imNum) = rot90(error,2);
    end
end

for filterNum = 1:numFilters
    for imNum = 1:numImages
        Wc_grad(:,:,filterNum) = Wc_grad(:,:,filterNum) + conv2(images(:,:,imNum),DeltaConv(:,:,filterNum,imNum),'valid');
    end
end
Wc_grad = (1./numImages) .* Wc_grad + lambda*Wc;

Step 2: Gradient Check

當時明明我的梯度下降沒法收斂，這一步居然通過了=。=

Step 3: Learn Parameters

這步比較簡單，根據UFLDL對隨機梯度下降的解釋，在minFuncSGD中加上衝量的影響就可以了：

%%% YOUR CODE HERE %%%
        velocity = mom*velocity+alpha*grad;
        theta = theta - velocity;

Step 4: Test

執行cnnTrain，最後準確率可以達到97%+

以上就可UFLDL上cnn的實現，最重要的是弄清楚每一層在每一個過程中需要進行的操作，我都總結在文章開頭的表格裡面了~matlab給我一個很大的感受就是矩陣的demension match，有時候知道公式是什麼樣的，但是實現起來要考慮矩陣的維度，兩個維度match的矩陣才能相乘或者相加，不過好處就是再不知道怎麼寫程式碼的時候可以結果維度match來寫程式碼。而且cnn debug起來真的好睏難，完全不知道是哪裡出了問題=。=

參考：