机器学习实战——支持向量机

目录

一、简介

1.概述

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

3.最大间隔

二、实验操作

1.简化版SMO算法

2.Platt SMO 算法

三、总结


一、简介

1.概述

支持向量机(SVM,也称为支持向量网络),支持向量机(support vector machines)是一种二分类模型,监督式学习 (Supervised Learning)的方法,主要用在统计分类 (Classification)问题和回归分析 (Regression)问题上。它的目的是寻找一个超平面来对样本进行分割,分割的原则是间隔最大化,最终转化为一个凸二次规划问题来求解。是机器学习中获得关注最多的算法没有之一。

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

机器学习实战——支持向量机

线性可分:可以很容易就在数据中给出一条直线将两组数据点分开
分隔超平面:将数据集分割开来的直线。数据点在二维平面上,分隔超平面就只是一条直线,但数据集是三维时,那么分隔超平面就是一个平面。依此类推,如果数据集是 N ( N ≥ 2 N(N\geq2 N(N≥2)维时,那么就需要一个 N − 1 N-1 N−1维的对象来分隔数据。该对象被称为超平面,也就是分类的决策边界。理想状态是分布在超平面一侧的所有数据都属于某个类别,而分布在另一侧的所有数据则属于另一个类别。
间隔:离分隔超平面最近的点,到分隔面的距离。间隔应该尽可能地大,这是因为如果我们犯错或者在有限数据上训练分类器的话,大的间隔可以增加分类器的鲁棒性。
支持向量:离分隔超平面最近的那些点。支持向量到分割面的距离应该最大化。
 

3.最大间隔

对一个数据点进行分类,当超平面离数据点的“间隔”越大,分类的确信度(confidence)也越大。所以,为了使得分类的确信度尽量高,需要让所选择的超平面能够最大化这个“间隔”值。这个间隔就是下图中的Gap的一半。

机器学习实战——支持向量机

 

分隔超平面的形式为:机器学习实战——支持向量机
点A到分隔超平面的法线长度:机器学习实战——支持向量机

最大化间隔的目标就是找出分类器定义中的w和b。为此,我们必须找到具有最小间隔的数据点,而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点,我们就需要对该间隔最大化。这就可以写作:

机器学习实战——支持向量机

直接求解上述问题相当困难,所以我们将它转换成为另一种更容易求解的形式。

机器学习实战——支持向量机 

约束条件为:

机器学习实战——支持向量机 

其中常数C 用于控制 “最大化间隔” 和 “保证大部分点的函数间隔小于1.0” 这两个目标的权重。在优化算法的实现代码中,常数C 是一个参数,因此可以通过调节该参数得到不同的结果。一旦求出了所有的 α,那么分隔超平面就可以通过这些 α 来表达。 

二、实验操作

1.简化版SMO算法

简化版SMO算法,省略了确定要优化的最佳 α 对的步骤,而是首先在数据集上进行遍历每一个 α ,再在剩余的数据集中找到另外一个 α ,构成要优化的 α 对,同时对其进行优化,这里的同时是要确保公式:机器学习实战——支持向量机

所以改变一个 α 显然会导致等式失效,所以这里需要同时改变两个 α 。 

代码如下:

from time import sleep
import matplotlib.pyplot as plt
import numpy as np
import random
import types

"""
函数说明:读取数据

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


"""
函数说明:随机选择alpha

Parameters:
    i - alpha
    m - alpha参数个数
Returns:
    j -
"""
def selectJrand(i, m):
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

"""
函数说明:修剪alpha

Parameters:
    aj - alpha值
    H - alpha上限
    L - alpha下限
Returns:
    aj - alpah值
"""
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

"""
函数说明:简化版SMO算法

Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
    无
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    #转换为numpy的mat存储
    dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
    #初始化b参数,统计dataMatrix的维度
    b = 0; m,n = np.shape(dataMatrix)
    #初始化alpha参数,设为0
    alphas = np.mat(np.zeros((m,1)))
    #初始化迭代次数
    iter_num = 0
    #最多迭代matIter次
    while (iter_num < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            #步骤1:计算误差Ei
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])
            #优化alpha,更设定一定的容错率。
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                #随机选择另一个与alpha_i成对优化的alpha_j
                j = selectJrand(i,m)
                #步骤1:计算误差Ej
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                #保存更新前的aplpha值,使用深拷贝
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                #步骤2:计算上下界L和H
                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
                #步骤3:计算eta
                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
                #步骤4:更新alpha_j
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                #步骤5:修剪alpha_j
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
                #步骤6:更新alpha_i
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #步骤7:更新b_1和b_2
                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
                #步骤8:根据b_1和b_2更新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("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
        #更新迭代次数
        if (alphaPairsChanged == 0): iter_num += 1
        else: iter_num = 0
        print("迭代次数: %d" % iter_num)
    return b,alphas

"""
函数说明:分类结果可视化

Parameters:
    dataMat - 数据矩阵
    w - 直线法向量
    b - 直线解决
Returns:
    无
"""
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

Parameters:
    dataMat - 数据矩阵
    labelMat - 数据标签
    alphas - alphas值
Returns:
    无
"""
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)

 结果为:

机器学习实战——支持向量机

 

2.Platt SMO 算法

在这两个版本(简化版和完整版)中,实现 α 的更改和代数运算的优化环节一模一样。在优化过程中,唯一的不同就是 选择 α 的方式。完整版的Platt SMO算法应用了一些能够提速的启发方法。
Platt SMO算法是通过一个 外循环 来选择第一个α值的,并且其选择过程会在两种方式之间进行交替:一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界 α中实现单遍扫描。而所谓非边界α指的就是那些不等于边界0或C的 α值。对整个数据集的扫描相当容易,而实现非边界 α 值的扫描时,首先需要建立这些 α值的列表,然后再对这个表进行遍历。同时,该步骤会跳过那些已知的不会改变的α值。
在选择第一个α值后,算法会通过一个 内循环 来选择第二个α值。在优化过程中,会通过 最大化步长 的方式来获得第二个α值。在简化版SMO算法中,我们会在选择j 之后计算错误率 E j。但在这里,我们会建立一个全局的缓存用于保存误差值,并从中选择使得步长或者说 E i − E j 最大的 α 值。
代码如下:

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算法
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算法
                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)

结果为:

机器学习实战——支持向量机 

 优化后的程序相较于之前用时更短。

三、总结

支持向量机有如下优缺点。

优点:

  1. 对于线性不可分的情况可以通过核函数,映射到高维特征空间实现线性可分。
  2. SVM学习问题可以表示为凸优化问题,因此可以利用已知的有效算法发现目标函数的全局最小值。而其他分类方法(如基于规则的分类器和人工神经网络)都采用一种基于贪心学习的策略来搜索假设空间,这种方法一般只能获得局部最优解。
  3. 小集群分类效果好。

缺点:

  1. SVM仅仅只限于一个二类分类问题,对于多分类问题解决效果并不好。
  2. 仅局限于小集群样本,对于观测样本太多时,效率较低。
  3. 寻求合适的核函数相对困难。

上一篇:数据可视化之MATPLOTLIB实战:PLT.POLAR()函数 绘制极线图 (转载)


下一篇:ML - 数据归一化