逻辑回归之基础知识及手写数字识别实例

作者：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 了，以 w 点乘 x 为自变量，函数图像如下:

(注：从图中可以看到 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 训练好了，千万不能忘记保存啊，这样，下次一个新数据过来的时候，我们就可以直接进行判断啦。