1. 支持向量机(Support Vector Machines)

SVM有很多实现,最流行的一种实现是序列最小优化(Sequential Minimal Optimization)算法
可以使用核函数(kernel)的方式将SVM扩展到更多数据集上
SVM的优缺点:

  • 优点:泛化错误率低,计算开销不大,结果易解释
  • 缺点:对参数调节和核函数的选择敏感,原始分类器不加修改仅适用于处理二类问题

2. 基于最大间隔分隔数据

可以用一条线将两类数据点分开,则这组数据线性可分,否则线性不可分
线性可分数据集中的分隔界限称为分隔超平面(separating hyperplane),也就是分类的决策边界,分布在超平面一侧的数据都属于某个类别,而另一侧的数据属于另一类别
我们希望找到离分隔超平面最近的点,确保它们离分隔面的距离尽可能的远,这里,点到分隔面的距离称之为间隔(margin)
支持向量(support vector)就是离分隔超平面最近的那些点

3. 寻找最大间隔

分割超平面的形式可以写成\(w^Tx+b\),要计算A点到分割超平面的距离,必须给出点到分隔面的法线长度,该值为\(|w^TA+b|/||w||\)

3.1 分类器求解的优化问题

输入数据给分类器会输出一个类别标签,使用类似Heaviside的函数,u<0时,f(u)输出-1,u>0为+1,即类别标签是-1和+1
当计算数据点到分隔面的距离并确定分隔面的放置位置时,间隔通过label \(w^Tx+b\)计算。现在的目标是找出分类器定义的w和b,为此需要找到具有最小间隔的数据点,而这些数据点就是前面提到的支持向量,一旦找到具有最小间隔的数据点,我们就需要对该间隔最大化,这可以写作:\(
\arg \max\limits_{w,b} \{ \min\limits_{n} (label\cdot (w^Tx+b)) \cdot \frac{1}{||w||}
\}\)
直接求解上式太难,需要转换形式。在上述优化问题中,给定了一些约束条件求最优值,约束条件就是label \(w^Tx+b\),因此该问题是一个带约束条件的优化问题,可以考虑拉格朗日乘子法。引入拉格朗日乘子,就可以基于约束条件来表示原问题,由于此处约束条件都是基于数据点的,可以将超平面写成数据点的形式:
$$ \max\limits_a \left[ \sum_{i=1}^m a-\frac{1}{2}\sum_{i,j=1}^m label^{(i)}\cdot label^{(i)}\cdot \alpha_{i} \cdot \alpha_{j}\langle x^{(i)},x^{(j)}\rangle
\right]$$
其约束条件为:\(\alpha \geq 0\)和\(\sum_{i=1}^m \alpha_i \cdot label^{(i)}=0\)
以上基于假设:数据必须100%线性可分,但实际并不是,这时引入”松弛变量(slack variable)”,来允许有些数据点可以处于分隔面的错误一侧,这样我们的优化目标就能保持不变,但是此时新约束规则为:
\(C\geq \alpha \geq 0,\)和\(\sum_{i=1}^m \alpha_i \cdot label^{(i)}=0\)
这里常数C用于控制”最大化间隔”和”保证大部分点的函数间隔小于1.0”这两个目标的权重,在优化算法的实现代码中将C作为一个参数,则可以通过调节这个参数得到不同的结果,一旦求出了所有的alpha,就可以用这些alpha里表示分隔超平面,SVM的主要工作就是求这些alpha

3.2 SVM应用的一般框架

  • 收集数据:可以使用任何方法
  • 准备数据:需要数值型数据
  • 分析数据:有助于可视化分隔超平面
  • 训练算法:SVM的大部分时间都源自训练,该过程主要实现两个参数的调优
  • 测试算法:十分简单的计算过程就可以实现
  • 使用算法:几乎所有的分类问题都可以使用SVM,值得一提的是,SVM本身是一个二类分类器,对多累问题应用SVM需要地代码做些修改

4. SMO高效优化算法

我们对上述最后的式子进行优化,其中一个是最小化的目标函数,一个是在优化过程中必须遵循的约束条件。以前,人们还在使用二次规划求解工具(quadratic solver)来求解上述优化问题,这种工具是一种用于在线性约束下优化具有多个变量的二次目标函数的软件。而这些二次规划求解工具则需要强大的计算能力的支撑,另外在实现上也十分复杂。所有需要做的围绕优化的事情就是训练分类器,一旦得到alpha最优值,我们就得到了分隔超平面并能够将其用于数据分类

4.1 Platt的SMO算法

SMO表示序列最小优化(Sequential Minimal Optimization)。Platt的SMO算法将大优化问题分解为多个小优化问题,这些小优化问题往往很容易求解,并且将它们进行顺序求解的结果与作为整体求解的结果是一致的,在结果完全相同的同时,SMO算法的求解问题短很多
SMO算法的目标是求出一系列alpha和b,一旦求出了这些alpha,就很容易计算出权重向量w并得到分隔超平面
SMO算法的工作原理是:每次循环中选择两个alpha进行优化处理,一旦找到一对合适的alpha,那么久增大其中一个同时减小另外一个,这里所谓的”合适”就是指两个alpha必须要符合一定的条件,条件之一就是这两个alpha必须在间隔边界之外,而第二则是这两个alpha还没有进行过区间化处理或者不在边界上

4.2 应用简化版SMO算法处理小规模数据集

Platt SMO算法的完整实现需要大量代码,于是从简化到完整逐渐递增
Platt SMO算法中的外循环确定要优化的最佳alpha对,而简化版却会跳过这一部分,首先在数据集上遍历每一个alpha,然后在剩下的alpha集合中随机选择另一个alpha,从而构建alpha对。这里有一点很重要,就是我们要同时改变两个alpha,之所以这样做是因为我们有一个约束条件:\(\sum \alpha_i \cdot label^{(i)}=0\)
由于改变一个alpha可能导致该约束条件失效,因此我们总是同时改变两个alpha
为此将构建一个辅助函数,用于在某个区间范围内随机选择一个整数。同时,我们也需要另一个辅助函数,用于在数值太大时对其进行调整,实现代码:
辅助函数:

def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

简化的SMO函数伪代码:

创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
    对数据集中的每个数据向量(内循环):
        如果该数据向量可以被优化:
            随机选择另外一个数据向量
            同时优化这两个向量
            如果两个向量都不能被优化,退出内循环
    如果所有向量都没被优化,增加迭代数目,继续下一次循环

简化版SMO算法:

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: print "L==H"; continue
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print "eta>=0"; continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas

4.3 利用完整Platt SMO算法加速优化

完整版Platt SMO的支持函数:

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag

def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek
        
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

完整Platt SMO算法中的优化例程

def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or\
       ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJ(i, oS, Ei)                                   
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H: print "L==H"; return 0
        eta = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - \
              oS.X[j,:]*oS.X[j,:].T
        if eta >= 0: print "eta>=0"; return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): 
             print "j not moving enough"; return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*\
                      (alphaJold - oS.alphas[j])
        updateEk(oS, i)                                      
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*\
             oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*\
             (oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*\
             oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*\
             (oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0

完整版Platt SMO的外循环代码:

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)
                print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print "iteration number: %d" % iter
    return oS.b,oS.alphas

5. 在复杂数据上应用核函数

核函数(kernel)将数据转换成易于分类器理解的形式

5.1 利用核函数将数据映射到高维空间

seventh
Figure 6-6中,数据点位于一个圆中,人类能够意识到这点,但对于分类器而言,他只能识别分类器的结果是大于0还是小于0,如果只在x、y轴构成的坐标系中插入直线进行分类的话,我们得不到理想结果。或许可以对圆中的数据进行某种形式的转换,从而得到新变量来表示数据。在此例中,我们将数据从一个特征空间转换到另一个特征空间,在新空间下我们很容易利用已有的工具对数据进行处理,数学家们喜欢将这个过程称为从一个特征空间到另一个特征空间的映射。通常情况下,这种映射会将低维特征空间映射到高维特征空间。
这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的,读者可以把核函数想象成一个包装器(wrapper)或接口(interface),它能把数据从一个很难处理的形式转换成另一个较容易处理的形式。换种简单但不确切的方式来理解,可以把它想象成另一种距离计算的方法,距离计算的方法有很多种,核函数一样具有多种类型。经过空间转换后,我们可以在高维空间中解决线性问题,也就等价于在低维空间中解决非线性问题。
SVM优化中的一个有点就是所有运算都可以写成内积(inner product)的形式,我们可以把内积运算转换为核函数,而不必简化处理。将内积替换为核函数的方式称为核技巧(kernel substation)

5.2 径向基核函数(radial basis function)

核函数不仅应用于SVM,还有其他很多机器学习算法也需要用到,下面来介绍一个流行的核函数:径向基核函数
径向基函数是一个采用向量作为自变量的函数,能够基于向量距离运算输出一个标量。这个距离可以是从<0,0>向量或其他向量开始计算的距离。径向基函数的高斯版本公式为:\(k(x,y)=\exp\left(\frac{-||x-y||^2}{2\sigma^2}\right)\)
\(\sigma\)是用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数
上述高斯核函数将数据从特征空间映射到更高维空间,具体说是映射到一个无穷维空间,使用高斯核函数会得到一个理想的结果。对于Figure 6-6中例子,可以发现直接度量数据点到原点距离即可,然而碰到一个非此类的数据集,我们便会陷入困境。在该数据集中,使用高斯核函数可以得到很好的结果,当然,该函数也可以用于许多其他的数据集,并得到较低的错误率。
添加函数kernelTrans(),对optStruct类进行修改:

def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

在线性核函数的情况下,内积计算在”所有数据集”和”数据集中的一行”这两个输入之间展开。在径向基核函数情况下,在for循环里对矩阵的每个元素计算高斯核函数的值,而在for循环结束后,我们将计算过程应用到整个向量中去。
为了使用该函数,前述innerL()和calcEk()需要做些修改,修改之后:

def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H: print "L==H"; return 0
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0: print "eta>=0"; return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0


def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

5.3 在测试中使用核函数

接下来我们构造一个对Figure 6-6中数据点进行有效分类的分类器,该分类器使用了径向基核函数。前面提到的径向基函数有一个用户定义的输入\(\sigma\),首先得确定它的大小,然后利用核函数构建分类器,整个函数如下:

def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) 
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd]                                    
    labelSV = labelMat[svInd];
    print "there are %d Support Vectors" % shape(sVs)[0]
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)

上述函数的可选输入参数即\(\sigma\)
可以尝试更换不同的k1参数以观察测试错误率、训练错误率、支持向量个数随k1的变化情况
支持向量的数目存在一个最优值,SVM的优点在于它能对数据进行高校分类。如果支持向量太少,就可能会得到一个很差的决策边界,如果支持向量太多,就相当于每次都利用整个数据集进行分类,这种分类方法称为k近邻
我们可以对SMO算法中的其他设置进行随意地修改或者建立新的核函数

6. 示例

手写识别问题, 使用SMO算法, 构建一个系统去测试手写数字上的分类器:

def img2vector(filename):
    returnVect = zeros((1,1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0,32*i+j] = int(lineStr[j])
    return returnVect

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)           #load the training set
    m = len(trainingFileList)
    trainingMat = zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     #take off .txt
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels    

def testDigits(kTup=('rbf', 10)):
    dataArr,labelArr = loadImages('trainingDigits')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] 
    labelSV = labelMat[svInd];
    print "there are %d Support Vectors" % shape(sVs)[0]
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadImages('testDigits')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)

基于SVM构建多类分类器已有很多研究和对比了,感兴趣可以查看C.W.Huset等人发表的一篇论文”A Comparison of Methods for Multiclass Support Vector Machines”

《机器学习实战》一书中,作者使用不同的σ值,并尝试了线性核函数,得到的结果是:
nineth
表中结果表明,径向基函数中的参数\(\sigma=10\)左右时,可以得到最小测试错误率,该参数值比前面例子中的取值1.3大得多,这些差距是由于数据不同造成的,手写识别的数据中有1024个特征,而这些特征的值可能高达1.0,而前面的例子中只有2个特征。如何取值呢?多次尝试可以得到比较好的结果,C的设置也会影响到分类的结果。
当然存在另外的SVM形式,把C同时考虑进了优化过程中,例如v-SVM,这里可以参考Sergios Theodoridis和Konstantinos Koutroumbas撰写的Pattern Recognition。
我们还会发现,最小的训练错误率并不对应于最小的支持向量数目,另一个值得注意的点是:线性核函数的效果并不是特别差,可以牺牲线性核函数的错误率来提高分类速度,