MachineLearning-SVM

发表于 | 分类于

Machine Learning-SVM

标签 ： MachineLearning 人工智能

算法概述

  支持向量机“SupporVectorMachines”属于有监督学习算法，是现有的最好的监督学习算法之一，是一种比较成熟的算法。其优势是泛化错误率较低，计算开销不大；但是其对于参数调节较为明显。在低维样本分界不规则时，对于核函数的选择也较为敏感。原始的SVM仅适用于处理二分类问题

算法思路

基础概念

  首先需要介绍的概念是“分割超平面”。对于样本集属性为二维的情况来说，分割超平面指的是二维坐标上的一条分割直线，用此直线将二维样本分为两类。对于更高维的情况来说，如果数据集属性为N维的，则需要一个N-1维的对象来将数据分割，此对象称为超平面，也就是分割界线。
  假如已经确定好分割超平面，对于样本来说，我们认为样本点距离分割超平面的距离越远，我们对于此预测的结果越有把握，而对于在此超平面附近的样本点我们则没有很大的把握预测结果是正确的。我们希望采用这种方法来构造分离器：样本集中样本点离决策边界越远，则结果越可信。由此，我们将问题变成寻找最大间隔。
  对于二分类的SVM来说，与Logistic回归不同，我们采用“1”和“-1”作为分类的标签，在这样的前提下 ，我们可以通过统一的公式来表示间隔或者数据点到分割平面的距离。

整体思路

  我们可以将分割超平面的形式表示为：
                  公式0
  在确定分割超平面表达式之后，由于我们要根据样本点到超平面的距离距离超平面的参数，我们需要定义样本点到分割超平面的距离。我们首先选择的距离成为函数距离，其公式为：
               公式1
  观察可以发现，此函数有以下优势：在分类正确的时候，函数距离为一个很大的正数，而在分类错误的时候。我们想找到分割超平面的一组参数样本集中里距离超平面距离最小的样本点离超平面距离最大。即使以下距离最大：
                  公式1
  这样就出现一个问题，如果我们只是将W变大的话，则该距离就会增大，但是这样的变化结果是所有样本点到超平面的距离都对应增大，显然是无效的，因此我们进行一步归一化处理，使用每个样本点的几何距离，公式如下：
               公式1
  进而最小化则为训练样本对于分界线的最小距离，为：
                   公式1
  在这样的处理之后，我们可以得到以下优化问题：
               公式1
  在经过一些处理之后，优化问题变为：
               公式1
  之后需要做的就是解此最优化问题。
  对于上述最优化问题，我们利用拉格朗日方法进行处理。将其转换为对偶优化问题，其表述如下：
             公式1
  SVM的工作主要就是求解参数alpha。

SMO算法

  目前我们常用的参数求最优化方法为SMO算法,SMO优化算法”SequentialMinimalOtimization”意为序列最小化优化。其主要思想是将最大优化问题分解为多个小优化问题来求解。小优化问题往往易于求解，并且它们求解的结果与将它们作为整体求解的结果完全一致。
  SMO算法的目标是求出一系列的alpha和b。其工作原理是：每次循环选择两个alpha进行优化处理。一旦找到一对合适的alpha，则增大其中一个同时减小另一个。“合适”的判定标准是两个alpha必须在间隔边界之外，并且没有进行过区间化处理或者不在边界上。即：
重复直到收敛 ：
  1. 选择某一对的ai和aj以在下次迭代中进行更新(这里需要选择那种能朝全局最大值方向最大程度靠近的一对值)
  2. 使用对应的ai和aj来重新优化W(a)，而保持其他的k值固定。

  其中参数alpha2迭代公式如下：
             公式2
  而alpha1参数更新由其与alpha2的约束条件可以求得：
                   公式1
  这两个alpha收敛之后，在继续其他alpha的迭代。只带到达最大迭代次数，算法结束。
  在训练完成进行预测时，利用如下转换即可：
               3

算法实现

  本文SVM实现的参考《MachineLearnigninAction》。
  首先实现简化版的SMO。其实现伪代码如下：
    创建一个alpha向量并将其初始化为0向量
    当迭代次数小于最大迭代次数时（外循环）：
      对数据集中的每个数据向量（内循环）：
      如果该数据向量可以被优化：
        随机选取另外一个数据向量
        同时优化这两个向量
        如果两个向量都不能被优化，退出循环
    如果所有向量都没有被优化，增加迭代数目，继续下一次循环
  实现代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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 ((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])
                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

  分类结果如图：
2
  然后是一个完整的SMO算法，是通过一个外循环来选择第一个alpha值的，并且选择过程在两种方法之间交替进行：一种是在所有数据集上进行单遍扫描，另一种是在非边界alpha中实现单遍扫描。在选择第一个alpha值会后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大步长的方法来获取第二个alpha值。
  完整SMO实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
class optStructK:
    def __init__(self,dataMatIn, classLabels, C, toler):  
        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)))

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 selectJK(i, oS, Ei):        
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   
            if k == i: continue 
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

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 = selectJK(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) #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                  
        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
    
def smoPK(dataMatIn, classLabels, C, toler, maxIter):    
    oS = optStructK(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    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 
        elif (alphaPairsChanged == 0): entireSet = True
        print "iteration number: %d" % iter
    return oS.b,oS.alphas

  结果如下：
1
  还有一个将alpha转化为w的函数实现如下：

1
2
3
4
5
6
7
def calcWs(alphas,dataArr,classLabels):
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

  具体计算的时候利用此函数得到参数w和b，带入样本属性，进行计算即可。

算法拓展

  与SVM紧密结合的是“核方法”，其作用是找到合适的核函数，将低维的数据投影到高维的空间内进行分类，用于处于在低维分类不明显的情况。算法对于核函数敏感。具体实现不再此做详细介绍。可以参考西瓜书和CS229的讲义，

算法小结

  支持向量机泛化错误率较低，具有良好的学习能力，且学习的结果有良好的推广性，在样本量比较少的时候，容易抓住数据和特征之间的非线性关系。