机器学习实战-第六章支持向量机

支持向量机（SVM）：

基本思想是求解能够正确划分训练数据集并且几何间隔最大的分离超平面

具体参考：

https://zhuanlan.zhihu.com/p/31886934 支持向量机的数学原理
https://zhuanlan.zhihu.com/p/26514613 KKT条件

SMO高效算法：
其中一个是对最小化的目标函数，一个是在优化过程中必须遵循的约束条件。

SMO算法的工作原理：
每次循环中选择两个Alpha进行优化处理，一旦找到一对合适的Alpha，那么就增大其中一个同时减小另一个，这里的合适是指符合两个条件
1.这两个Alpha必须要在间隔边界之外
2.这两个Alpha还没有进行过区间化处理或者不在边界上

在数据集上遍历每一个Alpha，然后在剩下的Alpha集合中随机选择另一个Alpha，从而构建Alpha对。

得满足条件∑αi*label = 0 所以需要同时改变两个Alpha

SMO算法中的辅助函数：

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

#i是第一个alpha的下标,m是所有alpha的数目
#只要函数值不等于输入值i,函数值就会随机选择
def selectJrand(i,m):
    j=i
    while(j==i):
        j = int(np.random.uniform(0,m))
    return j

#用于调整大于H或小于L的alpha值
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

#dataMatin数据集,classLabels类别标签,常数C,toler容错率,取消前最大的循环次数
def smoSimple(dataMatin, classLabels, C, toler, maxIter):
    dataMatrix = np.mat(dataMatin);labelMat = np.mat(classLabels).transpose()
    b = 0; m,n = np.shape(dataMatrix) #获取矩阵的行和列
    alphas = np.mat(np.zeros(m,1)) #alpha列矩阵,矩阵中元素都初始化为0
    iter = 0 #在没有任何alpha改变的情况下遍历数据集的次数
    #当变量达到maxIter时,函数结束运行并退出
    while (iter < maxIter):
        alphaPairsChanged = 0 #用于记录alpha是否已经进行优化
        for i in range(m):
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            #Ei是误差
            Ei = fXi - float(labelMat[i])
            #对Alpha值进行优化
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m) #选择Alpha
                fXj = float(np.multiply(alphas, labelMat).T,(dataMatrix*dataMatrix[j,:])) + b
                Ej = fXi - float(labelMat[j])
                #新旧值的比较
                alphaIold = alphas[i].copy() 
                alphaJold = alphas[j].copy()
                #保证alpha值不能等于0或者C
                #L和H用于将alpha[j]调整到0和C之间
                if(labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i]) #alpha小于0或者大于C被调整为0或C
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                #L等于H就不做任何调整
                if L==H:
                    print("L==H")
                    continue
                #eta是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
                #更新alpha_j
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                #步骤5：修剪alpha_j
                alphas[j] = clipAlpha(alphas[j],H,L)
                #检查alpha[j]是否有轻微改变,如果是就退出for循环
                if(abs(abs(alphas[j] - alphaJold) < 0.00001)):
                    print("j not moving enough")
                    continue
                #更新alpha_i
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #更新b1和b2   
                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[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T)
                #根据b1和b2更新b
                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))
        #在for循环之外,需要检查alpha值是否做了更新,如果有更新则将iter设为0后继续运行程序
        #只有在所有数据集上遍历maxIter次,且不再发生任何alpha修改之后,程序才会停止并退出while循环
        if(alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0 
            print("iterstaion number: %d" % iter)
    return b,alphas

def showClassifer(dataMat, w, b):
    #绘制样本点
    data_plus = []                                  #正样本
    data_minus = []                                 #负样本
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)              #转换为numpy矩阵
    data_minus_np = np.array(data_minus)            #转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)   #正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图
    #绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()

#计算w
def get_w(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
    w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()


if __name__ == '__main__':
    dataMat, labelMat = loadDataSet('testSet.txt')
    b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    w = get_w(dataMat, labelMat, alphas)
    showClassifer(dataMat, w, b)

加速优化的完整版Platt SMO算法：

在这两个版本（简化版和完整版）中，实现alpha 的更改和代数运算的优化环节一模一样。在优化过程中，唯一的不同就是 选择alpha 的方式。完整版的Platt SMO算法应用了一些能够提速的启发方法。

Platt SMO算法是通过一个 外循环 来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。

在选择第一个alpha值后，算法会通过一个 内循环 来选择第二个alpha值。在优化过程中，会通过 最大化步长 的方式来获得第二个alpha值。在简化版SMO算法中，我们会在选择j 之后计算错误率 Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说 Ei-Ej 最大的alpha 值。

import matplotlib.pyplot as plt
import numpy as np
import random

"""
   Parameters：
       dataMatIn - 数据矩阵
       classLabels - 数据标签
       C - 松弛变量
       toler - 容错率
   """
# 数据结构，维护所有需要操作的值（书上说是用于清理代码的数据结构）
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn#数据矩阵
        self.labelMat = classLabels#数据标签
        self.C = C#松弛变量
        self.tol = toler#容错率
        self.m = np.shape(dataMatIn)[0]#数据矩阵行数
        self.alphas = np.mat(np.zeros((self.m,1)))#根据矩阵行数初始化alpha参数为0   
        self.b = 0#初始化b参数为0
        #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。
        self.eCache = np.mat(np.zeros((self.m,2)))


"""
Parameters:
    fileName - 文件名
Returns:
    dataMat - 数据矩阵
    labelMat - 数据标签
"""
# 读取数据
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

"""
Parameters：
    oS - 数据结构
    k - 标号为k的数据
Returns:
    Ek - 标号为k的数据误差
"""
# 计算误差
def calcEk(oS, k):
    fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

"""
Parameters:
    i - alpha_i的索引值
    m - alpha参数个数
Returns:
    j - alpha_j的索引值
"""
# 函数说明:随机选择alpha_j的索引值
def selectJrand(i, m):
    j = i#选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

"""
Parameters：
    i - 标号为i的数据的索引值
    oS - 数据结构
    Ei - 标号为i的数据误差
Returns:
    j, maxK - 标号为j或maxK的数据的索引值
    Ej - 标号为j的数据误差
"""
# 内循环启发方式2
def selectJ(i, oS, Ei):
    maxK = -1; maxDeltaE = 0; Ej = 0#初始化
    oS.eCache[i] = [1,Ei]#根据Ei更新误差缓存
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]#返回误差不为0的数据的索引值
    if (len(validEcacheList)) > 1:#有不为0的误差
        for k in validEcacheList:#遍历,找到最大的Ek
            if k == i: continue#不计算i,浪费时间
            Ek = calcEk(oS, k)#计算Ek
            deltaE = abs(Ei - Ek)#计算|Ei-Ek|
            if (deltaE > maxDeltaE):#找到maxDeltaE
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej#返回maxK,Ej
    else:#没有不为0的误差
        j = selectJrand(i, oS.m)#随机选择alpha_j的索引值
        Ej = calcEk(oS, j)#计算Ej
    return j, Ej#j,Ej

"""
Parameters：
    oS - 数据结构
    k - 标号为k的数据的索引值
Returns:
    无
"""
# 计算Ek,并更新误差缓存
def updateEk(oS, k):
    Ek = calcEk(oS, k)#计算Ek
    oS.eCache[k] = [1,Ek]#更新误差缓存

"""
Parameters:
    aj - alpha_j的值
    H - alpha上限
    L - alpha下限
Returns:
    aj - 修剪后的alpah_j的值
"""
# 修剪alpha_j
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

"""
Parameters：
    i - 标号为i的数据的索引值
    oS - 数据结构
Returns:
    1 - 有任意一对alpha值发生变化
    0 - 没有任意一对alpha值发生变化或变化太小
"""
# 优化的SMO算法
# maxIter变量和函数smoSimple()中的作用有一点不同,后者当没有任何Alpha发生改变时会将整个集合的一次遍历过程计成一次迭代,而这里的一次迭代定义为一次循环过程,而不管该循环做了什么事
def innerL(i, oS):
    #步骤1：计算误差Ei
    Ei = calcEk(oS, i)
    #优化alpha,设定一定的容错率
    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)):
        #使用内循环启发方式2选择alpha_j,并计算Ej
        j,Ej = selectJ(i, oS, Ei)
        #保存更新前的aplpha值，使用深拷贝
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        #步骤2：计算上下界L和H
        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
        #步骤3：计算eta
        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
        #步骤4：更新alpha_j
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
        #步骤5：修剪alpha_j
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        #更新Ej至误差缓存
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("alpha_j变化太小")
            return 0
        #步骤6：更新alpha_i
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        #更新Ei至误差缓存
        updateEk(oS, i)
        #步骤7：更新b_1和b_2
        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
        #步骤8：根据b_1和b_2更新b
        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

"""
Parameters：
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
    oS.b - SMO算法计算的b
    oS.alphas - SMO算法计算的alphas
"""
# 完整的线性SMO算法
def smoP(dataMatIn, classLabels, C, toler, maxIter):
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)#初始化数据结构
    iter = 0#初始化当前迭代次数
    entireSet = True; alphaPairsChanged = 0
    #遍历整个数据集的alpha都没有更新或者超过最大迭代次数,则退出循环
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:#遍历整个数据集                           
            for i in range(oS.m):       
                alphaPairsChanged += innerL(i,oS)#使用优化的SMO算法,选择第二个Alpha
                print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:#遍历非边界值
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]#遍历不在边界0和C的alpha
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:#遍历一次后改为非边界遍历
            entireSet = False
        elif (alphaPairsChanged == 0):#如果alpha没有更新,计算全样本遍历
            entireSet = True 
        print("迭代次数: %d" % iter)
    return oS.b,oS.alphas#返回SMO算法计算的b和alphas

"""
Parameters:
    dataMat - 数据矩阵
    w - 直线法向量
    b - 直线解决
Returns:
    无
"""
# 分类结果可视化
def showClassifer(dataMat, classLabels, w, b):
    #绘制样本点
    data_plus = []#正样本
    data_minus = []#负样本
    for i in range(len(dataMat)):
        if classLabels[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)#转换为numpy矩阵
    data_minus_np = np.array(data_minus)#转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)#正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)#负样本散点图
    #绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()

"""
Parameters:
    dataArr - 数据矩阵
    classLabels - 数据标签
    alphas - alphas值
Returns:
    w - 计算得到的w
"""
# 计算w
def calcWs(alphas,dataArr,classLabels):
    X = np.mat(dataArr); labelMat = np.mat(classLabels).transpose()
    m,n = np.shape(X)
    w = np.zeros((n,1))
    for i in range(m):
        w += np.multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w


if __name__ == '__main__':
    dataArr, classLabels = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
    w = calcWs(alphas,dataArr, classLabels)
    showClassifer(dataArr, classLabels, w, b)

在复杂数据上应用核函数：

核函数的目的主要是为了解决非线性分类问题，通过核技巧将低维的非线性特征转化为高维的线性特征，从而可以通过线性模型来解决非线性的分类问题。

核函数：会将低维特征空间映射到高维空间。

SVM的优点在于它能对数据进行高效分类。如果支持向量太少，就可能会得到一个很差的决策边界；如果支持向量太多，也就相当于每次都利用整个数据集进行分类，这种分类方法称为k近邻。

核技巧：将内积替换成核函数的方法

核转换函数：

def kernelTrans(X, A, kTup):
    m,n = np.shape(X)
    K = np.mat(np.zeros((m,1)))
    if kTup[0] == 'lin': K = X * A.T#线性核函数,只进行内积。
    elif kTup[0] == 'rbf':#高斯核函数,根据高斯核函数公式进行计算
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = np.exp(K/(-1*kTup[1]**2))#计算高斯核K
    else: raise NameError('核函数无法识别')
    return K#返回计算的核K

基于SVM的手写识别问题：

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)           
    m = len(trainingFileList)
    trainingMat = np.zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     
        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   

"""
Parameters:
    kTup - 包含核函数信息的元组
Returns:
    无
"""
# 测试函数
def testDigits(kTup=('rbf', 10)):
    dataArr,labelArr = loadImages('trainingDigits')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10, kTup)
    datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
    svInd = np.nonzero(alphas.A>0)[0]
    sVs=datMat[svInd]
    labelSV = labelMat[svInd];
    print("支持向量个数:%d" % np.shape(sVs)[0])
    m,n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
    print("训练集错误率: %.2f%%" % (float(errorCount)/m))
    dataArr,labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
    m,n = np.shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1   
    print("测试集错误率: %.2f%%" % (float(errorCount)/m))


if __name__ == '__main__':
    testDigits()