2012-12-29

基于SVM的手写数字识别

前面两篇blog介绍了支持向量机SVM和SMO算法，这一篇就讲讲SVM的一个简单实际应用：使用svm实现一个简单的数字手写识别软件。首先要解决如何使用svm进行多类分类。

svm与多类分类¹

svm是一个二类的分类器，即它只回答属于正类还是负类的问题。而现实中要解决的问题，往往是多类的问题，比如文本分类，比如数字识别。如何由两类分类器得到多类分类器, 一般有三种方法：

一对多（一对其余）

比如我们有5个类别，第一次就把类别1的样本定为正样本，其余2，3，4，5的样本合起来定为负样本，这样得到一个两类分类器。通过类似的方法构造类别2、3、4、5的分类器，对于测试数据，每一个分类器都过一次，在那个分类器中，被判定为正，那么就认识它属于那个分类。但有两个问题还没有解决：一是被多个分类器判定为正，二是被所有分类器判定为负。第一种情况，可以简单地选择第一个被判定为正的分类，对二种情况，可以视为分类失败了。另外一个问题是，“其余”的那一类样本数总是要数倍于正类（因为它是除正类以外其他类别的样本之和嘛），这就人为的造成了的“数据集偏斜”问题。

一对一

这种方法需要构造n*(n-1)/2个分类器。简单地可以理解为所有分类1v1进行pk,具体做法，还是每次选一个类的样本作正类样本，而负类样本则变成只选一个类，(例如构造一个1 pk 2的分类器，它的正样本为1的分类，负样本为2的分类）这就避免了偏斜。因此过程就是算出这样一些分类器，第一个只回答“是第1类还是第2类”，第二个只回答“是第1类还是第3类”，第三个只回答“是第1类还是第4类”，如此下去，所以一共有5*(5-1)/2=10个分类器。

虽然分类器的数目多了，但是在训练阶段（也就是算出这些分类器的分类平面时）所用的总时间却比“一类对其余”方法少很多，在真正用来分类的时候，把一个测试数据扔给所有分类器，第一个分类器会投票说它是“1”或者“2”，第二个会说它是“1”或者“3”，让每一个都投上自己的一票，最后统计票数，如果类别“1”得票最多，就判这篇文章属于第1类。但这有个问题，分类器的数量是类别数量的平方,例如，类别数如果是1000，要调用的分类器数目会上升至约500,000个（类别数的平方量级）。

DAG SVM

这种方法是构造一个DAG SVM，（有向无环的svm)。还是像一对一方法那样来训练，只是在对一篇文章进行分类之前，先按照下面图的样子来组织分类器这样在分类时,我们就可以先问分类器“1对5”（意思是它能够回答“是第1类还是第5类”），如果它回答5，我们就往左走，再问“2对5”这个分类器，如果它还说是“5”，我们就继续往左走，这样一直问下去，就可以得到分类结果。

好处在哪？我们其实只调用了4个分类器（如果类别数是k，则只调用k-1个），分类速度飞快，且没有分类重叠和不可分类现象！缺点在哪？假如最一开始的分类器回答错误（明明是类别1的文章，它说成了5），那么后面的分类器是无论如何也无法纠正它的错误的（因为后面的分类器压根没有出现“1”这个类别标签），其实对下面每一层的分类器都存在这种错误向下累积的现象。不过不要被DAG方法的错误累积吓倒，错误累积在一对其余和一对一方法中也都存在，DAG方法好于它们的地方就在于，累积的上限，不管是大是小，总是有定论的，有理论证明。而一对其余和一对一方法中，尽管每一个两类分类器的泛化误差限是知道的，但是合起来做多类分类的时候，误差上界是多少，没人知道，这意味着准确率低到0也是有可能的，这多让人郁闷。

而且现在DAG方法根节点的选取（也就是如何选第一个参与分类的分类器），也有一些方法可以改善整体效果，我们总希望根节点少犯错误为好，因此参与第一次分类的两个类别，最好是差别特别特别大，大到以至于不太可能把他们分错；或者我们就总取在两类分类中正确率最高的那个分类器作根节点，或者我们让两类分类器在分类的时候，不光输出类别的标签，还输出一个类似“置信度”的东东，当它对自己的结果不太自信的时候，我们就不光按照它的输出走，把它旁边的那条路也走一走，等等。

高斯核函数

这个手写识别svm的核函数采用了高斯核函数。高斯核函数的公式如下

\[ K({x_i},{x_j})=e^{- \frac { {\left \| {x_i - x_j} \right \| }^2}{2\delta ^2}} \]

其中的径向基函数的宽度$\delta$对分类器的性能比较敏感，对取不同的$\delta$时，高斯核支持向量机的性能进行分析，若$0 $ ，则所有的训练样本点都是支持向量，且它们全部能被正确的分类，但容易出现“过学习”的现象，推广能力较差，对测试样本的错误识别率较高；若$\delta \to \infty$ ，高斯核支持向量机对所有样本一视同仁，推广能力或对测试样本的正确判别能力为零，即它把所有样本点判为同一类。实际上，当$\delta$取比训练样本点之间的平均距离小得多时，就能达到$0 $的效果；当$$取比训练样本点之间的平均距离大得多时，就能达到$$的效果。在确定高斯径向基函数的宽度$$时，最基本的方法是对$$取不同的值，然后分别采用支持向量机方法进行训练，选择最小分类错误率的一组$$参数。比较典型的方法有梯度下降法与交叉验证法。²

数据

我实现的这个数学手写识别的数据是“Machine Learning in Action”的第二章的数据，可以下载这本的source code获得。解压后的路径为machinelearninginaction/Ch02/digits.zip。

将digits.zip解压后就可以得到训练数据和测试数据。这些都是原始数据经过处理后的数据，是一个各个文件，一个文件是一个手写数据，文件名像这样“0_52.txt”，签名的0表示这个手写数据对应的数字，后面的数字52，是0这一类手写数据的文件编号。书写数据如下图所示，是一个32*32的黑白位图。

在处理手写数据时，例子中是将它们存放到一个32*32的矩阵中，也可以存放到一个有1024的一维数组。计算高斯核函数的${\| {x_i} - {x_j} \|} ^2 $如下，以一维数组为例：

\[ \begin{array}{l} {\left \| {x_i} - {x_j} \right \|}^2 = \sum\limits_{k=0}^{1023}{ {(x_i^k - x_j^k)}^2} \\ x_i=(x_i^0,x_i^1,\cdot \cdot \cdot,x_i^{1023}) \end{array} \]

实现

我实现了一个简单的基于svm的数字手写识别的python脚本，多类分类的模型采用的是一对多（一对其余）的方式，所以一共构造了10个分类器，分别对应0到9。算法采用的是前一篇blog讲的SMO，算法的实现也基本是一样的。下面是具体的代码，可以下载github上的代码和数据。

基于svm的数字手写识别view raw

#!/usr/bin/env python
# -*- coding: utf-8 -*-

### ###
# handwriting numbers OCR based on svm.
# data from "Machine Learning in Action" chapter 02.
# 基于SVM的手写数字的识别程序
# 数据：采用了《Machine Learning in Action》第二章的数据
### ###

import sys
import os
import math

class model:
    def __init__(self):
        self.a = []
        self.b = 0.0

class GV:
    def __init__(self):
        self.samples = []        # 样本数据
        self.tests = []          # 测试数据
        self.models = []         # 训练的模型
        self.diff_dict = []      # 用于缓存预测知与真实y之差Ei
        self.cur_mno = 0         # 当前正使用或训练的模型
        self.cache_kernel = []   # 缓存kernel函数的计算结果
        self.use_linear = False  # 是否使用线性核函数
        self.RBF_dlt = 10        # 径向基函数的宽度

    def init_models(self):
        for i in range(0, 10):
            m = model()
            for j in range(len(self.samples)):
                m.a.append(0)
            self.models.append(m)

    def init_cache_kernel(self):
        i = 0
        for mi in self.samples: 
            print i
            self.cache_kernel.append([])
            j = 0
            for mj in self.samples:
                if i > j:
                    self.cache_kernel[i].append(self.cache_kernel[j][i])
                else:
                    self.cache_kernel[i].append(kernel(mi,mj))
                j += 1
            i += 1

class image:
    def __init__(self):
        self.data = []
        self.num = 0
        self.label = []
        self.fn = ""
       
    def printself(self):
        print "data"
        for line in self.data:
            print line
        print "num", self.num
        print "label", self.label[gv.cur_mno]
        print "fn", self.fn 

# global variables
gv = GV()

def parse_image(path):
    img_map = []
    fp = open(path, "r") 
    for line in fp:
        line = line[:-2]
        img_map.append(line)
    return img_map

# load samples and tests
def loaddata(dirpath, col):
    files = os.listdir(dirpath)
    for file in files:
        img = image()
        img.data = parse_image(dirpath + file)
        img.num = int(file[0])
        img.fn = file
        col.append(img)

def kernel(mj, mi):
    if gv.use_linear == True:
        return kernel_linear(mj,mi)
    else:
        return kernel_RBF(mj,mi)

######
# 高斯核函数
######
def kernel_RBF(mj, mi):
    if mj == mi:
        return math.exp(0)
    dlt = gv.RBF_dlt
    ret = 0.0
    for i in range(len(mj.data)):
        for j in range(len(mj.data[i])):
            ret += math.pow(int(mj.data[i][j]) - int(mi.data[i][j]), 2)
    ret = math.exp(-ret/(2*dlt*dlt))
    return ret

######
# 线性
######
def kernel_linear(mj, mi):
    ret = 0.0
    for i in range(len(mj.data)):
        for j in range(len(mj.data[i])):
            ret += int(mj.data[i][j]) * int(mi.data[i][j])
    return ret


# g(x)
def predict(m):
    pred = 0.0
    for j in range(len(gv.samples)):
        if gv.models[gv.cur_mno].a[j] != 0:
            pred += gv.models[gv.cur_mno].a[j] * gv.samples[j].label[gv.cur_mno] * kernel(gv.samples[j],m)
    pred += gv.models[gv.cur_mno].b 
    return pred

# the same as predict(m), only with different parmaters
def predict_train(i):
    pred = 0.0
    for j in range(len(gv.samples)):
        if gv.models[gv.cur_mno].a[j] != 0:
            pred += gv.models[gv.cur_mno].a[j] * gv.samples[j].label[gv.cur_mno] * gv.cache_kernel[j][i]
    pred += gv.models[gv.cur_mno].b 
    return pred

# 决策函数对xi的预测值和真实值之差
def predict_diff_real(i):
    diff = predict_train(i)
    diff -= gv.samples[i].label[gv.cur_mno]
    return diff

# 优化计算Ei
def predict_diff_real_optimized(idx, i, new_ai, j, new_aj, new_b):
    diff = (new_ai - gv.models[gv.cur_mno].a[i])* gv.samples[i].label[gv.cur_mno] * gv.cache_kernel[i][idx]
    diff += (new_aj - gv.models[gv.cur_mno].a[j])* gv.samples[j].label[gv.cur_mno] * gv.cache_kernel[j][idx]
    diff += new_b - gv.models[gv.cur_mno].b
    diff += gv.diff_dict[idx]
    return diff


def init_predict_diff_real_dict():
    gv.diff_dict = []
    for i in range(len(gv.samples)):
        gv.diff_dict.append(predict_diff_real(i))

def update_diff_dict(i, new_ai, j, new_bj, new_b):
    for idx in range(len(gv.samples)):
        # 原来的函数
        # gv.diff_dict[idx] = predict_diff_real(idx)
        # 有优化后的
        gv.diff_dict[idx] = predict_diff_real_optimized(idx, i, new_ai, j, new_bj, new_b)

def update_samples_label(num):
    for img in gv.samples:
        if img.num == num:
            img.label.append(1)
        else:
            img.label.append(-1)

######
#  svmocr train
#  基于算法SMO
#  T: tolerance 误差容忍度(精度)
#  times: 迭代次数
#  C: 惩罚系数
#  Mno: 模型序号0到9
#  step: aj移动的最小步长
######
def SVM_SMO_train(T, times, C, Mno, step):
    time = 0
    gv.cur_mno = Mno
    update_samples_label(Mno)
    init_predict_diff_real_dict()
    updated = True
    while time < times and updated:
        updated = False
        time += 1
        for i in range(len(gv.samples)):
            ai = gv.models[gv.cur_mno].a[i]
            Ei = gv.diff_dict[i]

            # agaist the KKT
            if (gv.samples[i].label[gv.cur_mno] * Ei < -T and ai < C) or (gv.samples[i].label[gv.cur_mno] * Ei > T and ai > 0):
                for j in range(len(gv.samples)):
                    if j == i: continue
                    kii = gv.cache_kernel[i][i]
                    kjj = gv.cache_kernel[j][j]
                    kji = kij = gv.cache_kernel[i][j] 
                    eta = kii + kjj - 2 * kij 
                    if eta <= 0: continue
                    new_aj = gv.models[gv.cur_mno].a[j] + gv.samples[j].label[gv.cur_mno] * (gv.diff_dict[i] - gv.diff_dict[j]) / eta # f 7.106
                    L = 0.0
                    H = 0.0
                    a1_old = gv.models[gv.cur_mno].a[i]
                    a2_old = gv.models[gv.cur_mno].a[j]
                    if gv.samples[i].label[gv.cur_mno] == gv.samples[j].label[gv.cur_mno]:
                        L = max(0, a2_old + a1_old - C)
                        H = min(C, a2_old + a1_old)
                    else:
                        L = max(0, a2_old - a1_old)
                        H = min(C, C + a2_old - a1_old)
                    if new_aj > H:
                        new_aj = H
                    if new_aj < L:
                        new_aj = L
                    if abs(a2_old - new_aj) < step:
                        print "j = %d, is not moving enough" % j
                        continue

                    new_ai = a1_old + gv.samples[i].label[gv.cur_mno] * gv.samples[j].label[gv.cur_mno] * (a2_old - new_aj) # f 7.109 
                    new_b1 = gv.models[gv.cur_mno].b - gv.diff_dict[i] - gv.samples[i].label[gv.cur_mno] * kii * (new_ai - a1_old) - gv.samples[j].label[gv.cur_mno] * kji * (new_aj - a2_old) # f7.115
                    new_b2 = gv.models[gv.cur_mno].b - gv.diff_dict[j] - gv.samples[i].label[gv.cur_mno]*kji*(new_ai - a1_old) - gv.samples[j].label[gv.cur_mno]*kjj*(new_aj-a2_old)    # f7.116
                    if new_ai > 0 and new_ai < C: new_b = new_b1
                    elif new_aj > 0 and new_aj < C: new_b = new_b2
                    else: new_b = (new_b1 + new_b2) / 2.0

                    update_diff_dict(i, new_ai, j, new_aj, new_b)
                    gv.models[gv.cur_mno].a[i] = new_ai
                    gv.models[gv.cur_mno].a[j] = new_aj
                    gv.models[gv.cur_mno].b = new_b
                    updated = True
                    print "iterate: %d, changepair: i: %d, j:%d" %(time, i, j)
                    #break

# 测试数据
def test():
    recog = 0
    recog_correct = 0
    for img in gv.tests: 
        print "test for", img.fn
        for mno in range(10):
            gv.cur_mno = mno
            if predict(img) > 0:
                print mno
                print img.fn
                recog += 1
                if mno == int(img.fn[0]):
                    recog_correct += 1
                break
    print "recog:", recog
    print "recog_correct:", recog_correct
    print "total:", len(gv.tests)

def save_models():
    for i in range(10):
        fn = open("models/" + str(i) + "_a.model", "w")
        for ai in gv.models[i].a:
            fn.write(str(ai))
            fn.write('\n')
        fn.close()
        fn = open("models/" + str(i) + "_b.model", "w")
        fn.write(str(gv.models[i].b))
        fn.close()

def load_models():
    for i in range(10):
        fn = open("models/" + str(i) + "_a.model", "r")
        j = 0
        for line in fn:
            gv.models[i].a[j] = float(line)
            j += 1
        fn.close()
        fn = open("models/" + str(i) + "_b.model", "r")
        gv.models[i].b = float(fn.readline())
        fn.close()

if __name__ == "__main__":
    training = True
    loaddata("trainingDigits/", gv.samples)
    loaddata("testDigits/", gv.tests)    
    print len(gv.samples)
    print len(gv.tests)

    if training == True:
        gv.init_cache_kernel()
    gv.init_models()

    print "init_models done"

    T = 0.0001
    C = 10
    step = 0.0001
    gv.RBF_dlt = 8
    if training == True:
        for i in range(10):
            print "traning model no:", i
            SVM_SMO_train(T, 100, C, i, step)
        save_models()
    else:
        load_models()
        for i in range(10):
            update_samples_label(i)
    test()

训练数据一共是1934个，测试数据是946，一共是识别出来了914个，其中有911个识别正确。运行结果的截图如下：

在实现的过程中，第一个版本效率很低，主要是做了两个地方的改进，一是，缓存了核函数的运算结果，二是，计算预测值的与真实值之差$E\_i$，改用增量变化的方式，即每次只计算出$new\\_a\_i$, $new\\_a\_j$对$E\_i$的改变量。同样一开始，还有识别率低的问题，主要通过配置精度，惩罚系数，$a\_i$改变最小步长，这三个参数来提高识别率。但还是对高斯核函数的径向基函数的宽度$\delta$的调整最有效果，针对本例最优情况是当$\delta=10$左右。

不足

本文实现的例子，还有很多不足之处:

样本数据的偏移，训练时，正负样本的数量应该是相当的。本文中的例子，负样本是正样本的9倍。
高斯径向基函数的宽度，没有经过训练
选择第二个变量时，本文使用的是遍历方式，也可以改成寻找改变量最大的第二个参数，具体《统计学习方法》的P129
核函数计算时，效率比较低，主要影响测试分类时的速度，如果能用到位运算，应该效率会高很多。

其实上面的情况都还可以慢慢改进，最大的坑还是训练的时间太长了，我在虚拟机上装的ubuntu server运行的脚本，训练数据接近2000个，测试数据1000个，总共可能花了4个小时左右，时间太长。我之前的同事实现的ocr是居于online learing的，那么我们下一篇blog也来实现一个online learning的数字手写ocr，与今天这个对比一下。

主要参考了：SVM入门（十）将SVM用于多类分类。↩
张翔，肖小玲，徐光祐，一种确定高斯核模型参数的新方法，计算机工程，第33卷，第12期，2007年6月。↩

svm与多类分类1

一对多（一对其余）

一对一

DAG SVM

高斯核函数

数据

实现

不足

svm与多类分类¹