作者：daniel-D

原文：http://www.cnblogs.com/daniel-D/

  這學期 Pattern Recognition 課程的 project 之一是手寫數字識別，之二是做一個網站驗證碼的識別(鴨梨不小哇）。麵包要一口一口吃，先嚐試把模式識別的經典問題——手寫數字識別做出來吧。這系列部落格參考deep learning tutorial ，記錄下用以下三種方法的實現過程:

1. Logistic Regression - using Theano for something simple
2. Multilayer perceptron - introduction to layer
3. Deep Convolutional Network - a simplified version of LeNet5

  目的在於一方面從理論上幫自己理順思路，另一方面作為課程最後一課 presentation 的材料。這三種方法從易到難，準確率也由低到高，我們先從簡單的入手。

1. Binomial logistic regression model

  儘管線性分類器方法足夠簡單並且使用廣泛，但是線性模型對於輸出的 y 沒有界限，y 可以取任意大或者任意小（負數）的值，對於某些問題來說不夠 adequate, 比如我們想得到 0 到 1 之間的 probability 輸出，這時候就要用到比 linear regression 更加強大的 logistic regression 了。

y = w • x

  直覺上，一個線性模型的輸出值 y 越大，這個事件 P(Y=1|x) 發生的概率就越大。 另一方面，我們可以用事件的機率（odds）來表示事件發生與不發生的比值，假設發生的概率是 p ，那麼發生的機率（odds）是 p/(1-p) ， odds 的值域是 0 到正無窮，機率越大，發生的可能性越大。將我們的直覺與機率聯絡起來的就是下面這個（log odds）或者是 logit 函式 (有點牽強 - -!)：

  進而可以求出概率 p 關於 w 點乘 x 的表示：

(注：為什麼要寫出中間一步呢？看到第三部分的你就明白啦!)

  這就是傳說中的 sigmoid function

(注：從圖中可以看到 wx 越大，p 值越高，線上性分類器中，一般 wx = 0 是分介面，對應了 logistic regression 中 p = 0.5)

2. Parameter Estimation

  Logsitic regression 輸出的是分到每一類的概率，引數估計的方法自然就是最大似然估計 (MLE) 咯。對於訓練樣本來說，假設每個樣本是獨立的，輸出（標籤）為 y = {0, 1}，樣本的似然函式就是將所有訓練樣本 label 對應的輸出節點上的概率相乘, 令 p = P(Y=1|x) ,如果 y = 1, 概率就是 p， 如果 y = 0, 概率就是 1 - p ，（好吧，我是個羅嗦的傢伙), 將這兩種情況合二為一，得到似然函式：

  嗯？有連乘，用對數化為累加， balabala 一通算下來，就得到了對數似然函式為

  應用梯度下降法或者是擬牛頓法對 L(w) 求極大值，就可以得到 w 的估計值了。

3. Softmax regression

  這一小節旨在弄清 softmax regression 和 logistic regression 的聯絡，更多細節請參考 Andrew Ng 的英文資料，需要快速瀏覽也可以看看中文翻譯版本，或者 wiki 一下。

  logistic regression 在多類上的推廣又叫 softmax regression， 在數字手寫識別中，我們要識別的是十個類別，每次從輸入層輸入 28×28 個畫素，輸出層就可以得到本次輸入可能為 0, 1, 2… 的概率。得花點時間畫個簡易版本的，看起來更直觀:

  OK, 左邊是輸入層，輸入的 x 通過中間的黑線 w (包含了 bias 項)作用下，得到 w.x, 到達右邊節點， 右邊節點通過紅色的函式將這個值對映成一個概率，預測值就是輸入概率最大的節點，這裡可能的值是 {0, 1, 2}。在 softmax regression 中，輸入的樣本屬於第 j 類的概率可以寫成：

(注：對比第一部分提到過的中間一步，有什麼不同？)

注意到，這個迴歸的引數向量減去一個常數向量，會有什麼結果：

  沒有變化！這說明如果某一個向量是代價函式的極小值點，那麼這個向量在減去一個任意的常數向量也是極小值點，這是因為 softmax 模型被過度引數化了。(題外話：回想一下線上性模型中，同時將 w 和 b 擴大兩倍，模型的分界線沒有變化，但是模型的輸出可信度卻增大了兩倍，而在訓練迭代中， w 和 b 絕對值越來越大，所以 SVM 中就有了函式距離和幾何距離的概念)

  既然模型被過度引數化了，我們就事先確定一個引數，比如將 w1 替換成全零向量，將 w1.x = 0 帶入 binomial softmax regression ，得到了我們最開始的二項 logistic regression (可以動手算一算), 用圖就可以表示為

(注：虛線表示為 0 的權重，在第一張圖中沒有畫出來，可以看到 logistic regression 就是 softmax regression 的一種特殊情況)

  在實際應用中，為了使演算法實現更簡單清楚，往往保留所有引數，而不任意地將某一引數設定為 0。我們可以對代價函式做一個改動：加入權重衰減 (weight decay)。 權重衰減可以解決 softmax 迴歸的引數冗餘所帶來的數值問題。並且此時代價函式變成了嚴格的凸函式， Hessian矩陣變為可逆矩陣，保證有唯一的解。(感覺與線性分類器裡限制 ||w|| 或者設定某一個 w 為全零向量一樣起到了減參的作用，但是這個計算起來簡單清晰，可以用高斯分佈的 MAP 來推導，其結果是將 w 軟性地限制在超球空間，有一點 “soft” 的味道，個人理解^ ^)

  加入權重衰減後的代價函式是：

  等號右邊第一項是訓練樣本 label 對應的輸出節點上的概率的負對數，第二項是 weight decay ，可以使用梯度下降法和 L-BFGS 等演算法可以保證收斂到全域性最優解。總結起來，logistic regression 是 softmax regression 的一種特殊形式，前者是二類問題，後者是多類問題，前者的非線性函式的唯一確定的 sigmoid function (預設 label 為 0 的權重為全零向量的推導結果), 後者是 softmax, 有時候我們並不特意把它們區分開來。

好，基礎知識準備完畢，下面我們就要在數字手寫識別專案裡面實戰一下了。(^_^)

本文是應用 logsitic regression 模型對手寫數字識別的實現，整個程式是基於 MNIST 手寫數字資料庫進行 train， cross validate 和 test 的，如需下載 python 實現的原始碼，請點選這裡，你還可以在這裡下載資料集。 MNIST 資料庫由NYU 的 Yann LeCun 等人維護， Yann LeCun 自 1998 年以來就一直從事這方面的研究，實現的方法包括 linear classifier， K-NN， SVM， CNN 等，他提出的卷積神經網路是第一個真正多層結構學習演算法，利用空間相對關係減少引數數目以提高訓練效能，咱手寫數字識別的第三種方法就是基於這個演算法滴，Yan LeCun 同志曾經在深度機器學習大神 Geoffrey E. Hinton 底下做過博士後，現在也是 deep learning 的一個領軍人物了。

  言歸正轉，要用 python 實現這個專案還得用到 python 裡面一個比較特殊的 deep learning 的庫——Theano, 初次接觸這個庫，理解起來還需要一點時間，比如說 GPU 加速處理時，你需要將向量塊結構的變數轉換為 shared variables, 比如說類似於函式作用的 Graph structure, 學習曲線稍顯陡峭，如果你現在程式設計暫時用不到這些，直接跳過也行。

1. MNIST 相關

  動手之前，我們要先了解一下 MNIST 的相關情況。MNIST 包含了 7w 張 28×28 pixels 大小、數字 size 已經歸一化 (標準是什麼？是每類數字外面輪廓框框的大小嗎？)、中心化( maybe 框框的中心)的圖片。 deep learning tutorial 已經把它分為 train set, validation set and test set 3 個部分，分別包含了 5w、 1w 和 1w 張圖片。你如果一定要自己親眼瞅一瞅，跑一跑下面的程式碼，看一看 train set 上第 0 張圖片長啥樣：

 1 import cPickle 2 import gzip 3  4 import numpy 5 import Image 6  7 f = gzip.open('../data/mnist.pkl.gz', 'rb') 8 train_set, valid_set, test_set = cPickle.load(f) 9 f.close()10 11 print "The 0'th label:", train_set[1][0]12 lst = numpy.asarray(train_set[0][0], dtype=numpy.float32)13 img = Image.fromstring(mode='F', size=(28, 28), data=lst*255)14 img.show()

  不出意外的話，第 0 張訓練圖片應該是 5, label 也是 5:

The 0'th label: 5

2. 梯度下降

  首先，本次實驗的 logistic regression 框架可以用下面這張圖來表示：

x: 28*28 個 pixels作為 784 個輸入節點 w: 輸入節點的權重 b: bias 項的權重 0~9: 輸出節點

此次實驗的優化方法採用的是梯度下降法。梯度下降法中我們熟知的是標準的梯度下降演算法 (ordinary gradient descent)， 即每次迭代需要計算基於所有 datapoint 的 loss function 的 gradient， 計算量很大，速度較慢，而且因為標準誤差曲面只有一個，如果存在多個極小值點，此演算法容易陷於區域性極小值。同時還有另外一種比較高效的方法是 SGD (stochastic gradient descent)， 這種方法 estimating the gradient from just a few examples at a time instead of the entire training set. 由於每次計算 gradient 的時候每個誤差曲面是不同的，它們不一定具有相同的極小值，所以有時可以避免區域性極小值 （這是我的直觀理解）。A stationary point with respect to the error funciton for the whole data set will generally not be a stationary point for each data point individually (Bishop PRML). 在這個程式中，考慮到 shared variable 的特點，我們採用的是類似 SGD 的 Minibatch SGD ，每次迭代是基於好幾百個 examples 的梯度計算。

  同時，由於我們採用的是梯度下降策略，並不要求 Hessian 矩陣是可逆的，所以損失函式無需 weight decay 項，直接最小化最大似然函式的負對數即可，當然，這種方法可能導致的後果是最優解不是唯一的，並且容易產生 overfitting （這裡採用了 early-stopping 方法解決這個問題）。對於某一次 batch of train set 來說， Loss Function 可以寫成：

  Logistic Regression 模型建立的程式碼如下：

 1 #################### 2 #Build Actual Model# 3 #################### 4  5 # allocate symbolic variables for the data 6 index = T.lscalar()  # index to a [mini]batch 7 x = T.matrix('x')  # the data is presented as rasterized images 8 y = T.ivector('y')  # the labels are presented as 1D vector of [int] labels 9 10 # 建立之前定義的 LR 類，這裡沒有寫出來。11 # 類裡面初始化了 w， b 為全零向量12 # 按照 LR 的結構建立了從 inputs 到 outputs 的 graph structure。13 classifier = LogisticRegression(input=x, n_in=28 * 28, n_out=10)14 15 # 輸入 batch 的 index， 函式按照 index 將資料塊傳遞給 x， 16 # x 作為 classfier 的預設引數進入 LogisticRegression 類處理流程17 # 最後輸出錯誤分類的點的數目18 # 整個過程就像 pipes19 test_model = theano.function(inputs=[index],20         outputs=classifier.errors(y),21         givens={22             x: test_set_x[index * batch_size: (index + 1) * batch_size],23             y: test_set_y[index * batch_size: (index + 1) * batch_size]})24 25 validate_model = theano.function(inputs=[index],26         outputs=classifier.errors(y),27         givens={28             x: valid_set_x[index * batch_size:(index + 1) * batch_size],29             y: valid_set_y[index * batch_size:(index + 1) * batch_size]})30 31 # 求梯度，直接呼叫函式32 g_W = T.grad(cost=cost, wrt=classifier.W)33 g_b = T.grad(cost=cost, wrt=classifier.b)34 35 updates = [(classifier.W, classifier.W - learning_rate * g_W),36            (classifier.b, classifier.b - learning_rate * g_b)]37 38 cost = classifier.negative_log_likelihood(y)39 40 # 資料先利用 givens 得到資料塊41 # 然後進入類求得 cost 函式42 # 最後利用 update 將權重更新43 train_model = theano.function(inputs=[index],44         outputs=cost,45         updates=updates,46         givens={47             x: train_set_x[index * batch_size:(index + 1) * batch_size],48             y: train_set_y[index * batch_size:(index + 1) * batch_size]})

3. Early Stopping

Early stopping是應對 overfitting 的方法之一。主要思路如下：先用 train set 將 model 進行訓練，然後用得到的 model 來預測 validation set， 預測效果可以用 error 來表示。如果 error 在減小，說明我們的模型還可以繼續訓練，當 error 增大的時候，很有可能我們的 model 就 overfitting 了，這時候優化演算法就應該 halts 了。但是有一個問題是 “increasing validation error” is ambiguous， 很有可能是整體先下降再上升，但是在區域性上表現為 up and down， 就像股票走勢一樣。要解決這個問題，可以一開始我們就把 model 訓練到全域性極小值，然後對這個過程中的每一次迭代進行 validation error 的計算，這樣做雖然很安全，但是損失了 early stopping 快速的優點。這裡，我們採用的方法是：計算每一次迭代中 validation error， 如果比上次至少降低了 0.05%， 說明效果可以，就可以將迭代的限制次數增加到本次迭代次數的兩倍。顯然，這是一個基於經驗的辦法。

  訓練模型的重要引數：

1. epoch: 預設在整個 train set 層次上訓練輪數低於 1000
2. practice: 預設在 batch 層次上最少的迭代次數是 5000, 如果 model 在 validation set 上表現效果好，可以增加到目前迭代次數的 2 倍
3. validation_frequency: 每隔多少個 batch 評價一下 model， model 效果提升了，就可以增加 batch 的迭代限制次數。

  最後，由於 model 是 基於 validation set 上的最小錯誤率選出來的，因此，validation set 對於這個 model 來說是有偏的。換句話說， 選出來的模型會比較契合你現在的 validation set，所以不能用它來表示你模型的準確度，於是乎，我們又劃分出一個 model 從來沒有見過的 test set，用它的 error 來表示最終的預測能力。根據 Andrew Ng 的課程， train set : validation set : test set = 6 : 2 : 2, 這裡，我們用的是 5 : 1 : 1.

  模型訓練的程式碼如下:

1 #############2 #Train Model#3 #############4 epoch = 05 while (epoch < n_epochs) and (not done_looping):6         epoch = epoch + 17         for minibatch_index in xrange(n_train_batches):8                 minibatch_avg_cost = train_model(minibatch_index)9 ....

 由於程式碼比較簡單，這裡不做詳細解析了，詳細步驟請參考這裡的最後。

  用 logistic regression 的方法，可以將 test set 上的錯誤率降低到 7.489583 %。

4. Load and Save Models

  辛辛苦苦把 model 訓練好了，千萬不能忘記儲存啊，這樣，下次一個新資料過來的時候，我們就可以直接進行判斷啦。