支持向量机

机器学习实战

支持向量机（SVM）被认为是最好的现成分类器，意味着在数据上应用基本形式的SVM分类器就可以得到低错误率的结果。SVM能够对训练集之外的数据点做出很好的分类决策。这里介绍最流行的一种是实现即序列最小化算法。

基于最大间隔分隔数据

A
B

上面两组图中能够都画出一条直线将圆形点和方形点进行区分？
B组图中能够很容易的画出这样的一条直线将数据集分开，这样的数据叫做线性可分数据。这样的一条直线称为分隔超平面。在上线给出的例子当中，由于数据点是二维的，所以分隔超平面是一条直线。如果是数据是三维的则这分隔超平面是一个面。我们希望找到距离分隔超平面最近的点，确保它们离分隔面的距离尽可能给你的远。这里点到分隔面的距离被称为间隔。

支持向量机 就是离分隔超平面最近的那些点。

我们就是要找到最大化支持向量到分隔面的距离，需要找到此问题的优化求解方法。

寻找最大间隔

图3

图3中和都是可以分隔数据的超平面，但是怎么求解最佳分隔直线呢？
分隔超平面的形式可以定义为的形式，要计算图中点到超平面的距离需要给出点到分隔面的法线或垂线的长度，即。
其中

W:
X: 训练实例
b: bias

图4

假设现在是二维的情况下
X = (x1, X2)。
则超平面方程是 ：

则所有超平面右上方的点满足：
所有超平面左下方的点满足：

这样我们可以定义个函数是的,当u<0时f(u)输出-1，反之输出+1。当计算数据点到分隔面的距离并确定分隔面的位置时，间隔通过来计算。如果数据点处于正方向并且里分隔超平面很远的位置时，会是一个很大的正数，同时也会是一个很大的正数。而如果数据点处于负方向并且离分隔超平面很远的位置时，此时由于类别标签为-1，怎任然是个很大的正数。
现在需要找出其中w和b。为此我们需要找到具有最小间隔的数据点，而这些数据点也是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。

• 代码示例

//引入numpy
from numpy import *
//读取文件数据
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;
    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



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

'''
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.在对整个集合顺序遍历。用于记录alpha是否优化
        alphaPairsChanged = 0
        for i in range(m):
            #计算fXi这就是我们预测的类别 然后基于这个实例的预测结果和真是结果的对比 就可以计算Ei误差
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
            #如果误差很大就需要对该实例对应的alpha值进行优化
            #if语句中不管是正间隔还是负间隔都会被测试，并且在该if语句中，也同时检查alpha的值，以保证其不能等于0或C。
            #由于后面alpha小于0或者大于C时将会被调整为0或C，所有一旦在该if语句中等于这两个值的话，就不需要再优化了
            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])
                #计算两个alpha的误差
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                #计算L和H 用于将alpha[j]调整到0到C之间
                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
                #alpha[j]的最优修改量
                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] 如果轻微改变则结束循环
                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

if __name__ == '__main__':
    dataMat,labelMat = loadDataSet('svn.txt')
    b,alphas = smoSimple(dataMat,labelMat,0.6,0.001,40)
    for i in range(100):
        if alphas[i]>0:
            print(dataMat[i],labelMat[i])

[4.658191, 3.507396] -1.0
[3.457096, -0.082216] -1.0
[6.080573, 0.418886] 1.0

• sklearn示例

from sklearn import svm
if __name__ == '__main__':
    dataMat,labelMat = loadDataSet('svn.txt')
    #SVC有很多参数提供调试选择，达到最优的解
    clf = svm.SVC(kernel='linear')
    clf.fit(dataMat,labelMat)
    print(clf)
    print(clf.support_vectors_)

SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma='auto', kernel='linear',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
[[ 4.658191  3.507396]
 [ 3.457096 -0.082216]
 [ 6.080573  0.418886]]

利用完整Platte SMO算法加速优化

'''
用一个对象来保存一些重要的值
'''
class optStruct:
    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)))
'''
对于给定的alpha值 calcEk能够计算E值并返回 真实值与预测值之间的误差
'''
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
'''
用于选择第二个alpha值以保证在每次优化中采用最大的步长
'''
def selectJ(i,oS,Ei):
    maxK = -1;maxDeltaE = 0;Ej = 0
    oS.eCache[i] = [1,Ei]
    #返回非零E值对应的alpha值
    validEcacheList = nonzero(oS.eCache[:,0].A)[0] #.A是什么操作
    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 = 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
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):  # full Platt SMO
    oS = optStruct(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  # toggle entire set loop
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas