机器学习第六章笔记——支持向量机

最新推荐文章于 2023-05-07 20:54:38 发布

jgq1466693

最新推荐文章于 2023-05-07 20:54:38 发布

阅读量581

点赞数

文章标签：机器学习支持向量机 jvm

本文链接：https://blog.csdn.net/jgq1466693/article/details/126902395

版权

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

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

一、基于最大间隔分隔数据

支持向量机的优缺点：

优点：泛化错误率降低，计算开销不大，结果易解释

缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于处理二类问题

适用数据类型：数值型核标称型数据

观察图1和图2重点数据集：

通过上面观察上面两组数据集，可以看出对于线性可分数据集，可以轻易的画出一条直线将其分开，这条直线又叫分隔超平面，因为给出的例子是二维的，所以此时的超平面就是一条直线，若数据是三维的，那么此时用来分割数据的就是一个平面了。将N-1维的对象称为超平面，也就是分类的决策边界，分布在超平面一侧的所有数据都属于某种类别，另一侧的所有数据属于另一个类别。

如果数据点离决策边缘越远，那么最后的预测结果也就也可信。因此希望找到离分隔超平面最近的点，确保它们离分隔面的距离尽可能远，这里点到分隔面的距离被称为间隔。同时也希望间隔尽可能的大，这样即使犯错或者在有限数据上训练分类器时，也希望分类器能够足够的健壮。

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

二、寻找最大间隔

观察图3，分隔超平面的形式可以写成 $w^{T}x+b$ ，要计算点A到分隔超平面的距离，就必须要给出点到分隔面的法线或垂线的长度，该值为： $\left | w^{T}A+b \right | / \left \| w \right \|$ ，常数b类似于Logistic回归中的截距 $w_{0}$ ，这里的向量w和常数b一起描述了所有给的数据的分隔线或超平面。

2.1 分类器求解的优化问题

输入数据给分类器会输出一个类别标签，相当于一个类似于Sigmoid的函数在作用。利用单位阶跃函数对 $w^{T}x+b$ 作用，得到对应的函数值，以 $w^{T}x+b$ 为u，u<0时，输出-1，反之输出+1。至于为何选择-1和+1，是因为-1和+1仅仅相差一个符号，方便数学上的处理。通过一个统一的公式来表示间隔或数据点到分隔超平面的距离。

计算数据点到分隔面的距离并确定分隔面的放置位置时，间隔通过label*( $w^{T}x+b$ )来计算，如果数据点处于正方向并且离分隔超平面很远的位置时， $w^{T}x+b$ 将会是一个很大的正数。当数据点位于负方向并且离分割超平面很远的位置时，由于类别标签为-1， label*( $w^{T}x+b$ )仍会是一个正数。

如何找出分类器中定义的w和b是当前目标。找到具有最小间隔的数据点，这些数据点就是前面提到的支持向量。对间隔最大化，可以将其写做：

对于上式直接求解很困难，因此将其转换成另一种更容易求解的形式，观察上式中大括号内的部分，固定其中一个因子并最大化其他因子。令所有支持向量的 label*( $w^{T}x+b$ )都为1，只有那些离分隔超平面最近的点得到的值才为1.离超平面越远的数据点， label*( $w^{T}x+b$ )越大。

通过引入拉格朗日乘子，可以基于约束条件来表达上面的公式，这里的约束条件是

label*( $w^{T}x+b$ )>=1.0

约束条件是基于数据点的，因此可以将超平面写成数据点的形式，优化目标函数可以写为：

约束条件：

但是这都是基于数据必须100%线性可分的假设条件下进行。但就目前而言，几乎所有的数据都做不到这一点。引入松弛变量，允许有些数据点可以处于分隔面的错误一侧，这样优化目标就能保持仍然不变此时的新约束条件就为：

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

2.2 SVM应用的一般框架

一般流程：

1、收集数据：可以适用任意方法

2、准备数据：需要数值型数据

3、分析数据：有助于可视化分隔超平面

4、训练算法：SVM的大部分时间都源自训练，该过程主要实现两个参数的调优

5、测试算法：十分简单的计算过程就能够实现

6、适用算法：几乎所有分类问题都可以使用SVM，其本身是一个二类分类器，多多类问题应用SVM需要对代码进行修改。

三、SMO高效优化算法

3.1 Platt的SMO算法

SMO表示序列最小优化，该算法是将大优化问题分解成很多小优化问题进行求解，这些小优化问题往往很容易进行求解，并且对它们进行顺序求解的结果于将它们作为整体来求解的结果完全一致。但在结果完全一致的同时，SMO算法所需时间短的多。工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha，就增大其中一个同时减小另一个。所谓合适就是指两个alpha必须要符合在间隔边界之外以及两个alpha还没有进行区间化处理或者不在边界上。

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

简化版的SMO算法较完整的SMO算法而言，虽然代码减少了很多但其执行速度很慢。首先要在数据集上遍历每一个alpha，然后在剩下的alpha集合中选择另一个alpha，从而构建alpha对。

约束条件：

因为改变一个alpha可能会导致该约束条件失效，因此总是同时改变两个alpha。

构建一个辅助函数，用于在某个区间范围内随机选择一个整数，同时，需要另一个辅助函数，用于在数值太大时对其进行调整，下面的程序清单给出了这两个函数的实现。创建svmMLiA.py。

SMO算法的前期准备：

def loadDataSet(fliename):
    dataMat = []; labelMat = []
    fr = open(fliename)
    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))  # 返回0-m之间的随机浮点数
    return j

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

loaddataset()函数得到文本文件的每行的类标签以及整个数据矩阵。函数selectJrand()有两个参数值，其中i是第一个alpha的下标，m是所有alpha的数目，只要输入函数值不为输入值i，函数就会随机选择。

运行结果如下，采用的标签是-1，+1。

第一个版本的SMO算法：

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)  # 吗，n分别为dataMatrix的行和列
    # 初始化alpha参数，设为0
    alphas = np.mat(np.zeros((m, 1)))  # 生成m行1列的0矩阵
    # 初始化迭代次数
    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)):
                # toler容错率，c常数
                # 随机选择另一个与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

在每次的循环中，将alphaPairsChanged设置为0，对整个集合顺序遍历。变量alphaPairsChanged用以记录alpha是否已经完成了优化。fxi是预测的类别，基于这个预测的结果和真实结果的对比，计算误差Ei。当误差很大时，可以对该数据示例所对应的alpha值进行优化。

在if语句中不管是正间隔还是负间隔都会被测试。同时还要检测alpha值，以保证其不能等于0或c。之后就是选择第二个alpha值，即代码中的alpha[j]和第一个alpha值的操作类似，计算误差，可以通过copy()的方法实现，将老的alpha值与新的alpha值进行比较。计算L和H，用于调整alpha，使其在0-C之间。当L与H相等时，不做任何改变。Eta是alpha[j]的最优修改量，当eta为0，退出for循环的迭代过程。这里就能够得到新的alpha[j]，利用L、H对其进行调整。对alpha[i]和alpha[j]进行调整后给这两个alpha值设置一个常数项b。

测试效果：

dataArr, labelArr = svmMLiA.loadDataSet('D:\\learning\\testSet.txt')
b, alphas = svmMLiA.smoSimple(dataArr, labelArr, 0.6, 0.001, 40)

完整代码：

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


def loadDataSet(fliename):
    dataMat = []; labelMat = []
    fr = open(fliename)
    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))  # 返回0-m之间的随机浮点数
    return j

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


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)  # 吗，n分别为dataMatrix的行和列
    # 初始化alpha参数，设为0
    alphas = np.mat(np.zeros((m, 1)))  # 生成m行1列的0矩阵
    # 初始化迭代次数
    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)):
                # toler容错率，c常数
                # 随机选择另一个与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



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()


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('D:\\learning\\testSet.txt')
    b, alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    w = get_w(dataMat, labelMat, alphas)
    showClassifer(dataMat, w, b)

最终运行结果：

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

简化版SMO算法处理小规模数据集并没有太大的问题，但在更大的数据集上的运行速度就会变得很慢。完整版的SMO算法应用了一些能够提速的启发方法。

Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替，一中是在所有数据集上进行单遍扫描，另一种是在非边界alpha中实现单遍扫描。非边界alpha指那些不等于0或C的alpha值。优化过程中，通过最大化步长的方式获得第二个alpha值。通过建立一个全局的缓存用于保存误差值，并从中选择使得步长或者Ei-Ej最大的alpha值。

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


# 数据结构，维护所有需要操作的值（书上说是用于清理代码的数据结构）
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)))


# 读取数据
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 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


# 函数说明:随机选择alpha_j的索引值
def selectJrand(i, m):
    j = i  # 选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return 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


# 计算Ek,并更新误差缓存
def updateEk(oS, k):
    Ek = calcEk(oS, k)  # 计算Ek
    oS.eCache[k] = [1, Ek]  # 更新误差缓存


# 修剪alpha_j
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# 优化的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


# 完整的线性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


# 分类结果可视化
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()


# 计算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('D:\\learning\\testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
    w = calcWs(alphas, dataArr, classLabels)
    showClassifer(dataArr, classLabels, w, b)

最终运行结果为：

结果图如下：

五、在复杂数据上应用核函数

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

对于分类器而言，它只能够识别其结果是大于0还是小于0。若只在x和y轴构成的坐标系中插入直线进行分类，就并不会得到理想的结果。考虑对图5中的圆中数据进行某种形式的转换，从而得到某些新的变量来表示数据。这种情况下更容易得到大于0或小于0的测试结果。这里实现了将数据从一个特征空间转换到另一个特征空间，也可以将这个过程称为从一个特征空间到另一个特征空间的映射。

这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。可以将核函数理解为一个包装器或者是接口。可以把数据从某个很难处理的形式转换成另一个较容易处理的形式。

SVM优化中一个特别好的地方就是，所有的运算都可以写成内积的形式。将内积替换成核函数，而不必做简化处理。将内积替换成核函数的方式就称为核技巧。或者核”变电“。

5.2 径向基核函数

径向核函数是SVM中常用的一个核函数，径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。径向基函数的高斯版本，其具体公式为：

$\sigma$ 是用户定义的用于确定到达率或者说函数值跌落到0的速度参数。

上述高斯核函数将数据从其特征空间映射到更高维的空间，具体来说就是映射到一个无穷维的空间。在svmMLiA.py文件中添加一个函数并稍微的做一些修改。

def kernelTrans(X, A, kTup):
    m,n = shape(X)
    K = max(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 = np.exp(K / (-1 * kTup[1] ** 2))
    else:
        raise NameError('meet problem')
    return K


class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        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, 2)))
        self.b = 0
        self.ecache = np.mat(np.zeros((self.m, 2)))
        self.K = np.mat(np.zeros(self.m, self.m))
        for i in range(self.m):
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)

kTup是一个包含核函数信息的元组，在初始化方法结束时，矩阵K先被构建，再通过调用函数kernelTrans()进行填充。计算矩阵K时，多次调用了kernelTrans()，函数有3个输入参数：2个数值型变量和一个元组。元组kTiup给出的是核函数的信息。元组的第一个参数是描述所用核函数类型的一个字符串，其他参数是核函数可能需要的可选参数。在线性核函数的情况下，内积计算在”所有数据集“和数据集中的一行这两个输入之间展开，径向基函数的情况下，在for循环中对于矩阵的每一个元素计算高斯函数的值，当for循环结束后，将计算过程应用到整个向量上去。对于无法识别的元组，程序会抛出异常。

为了能够使用核函数，对lnnerL()和calcEk()进行修改:

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

# 优化的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, j] - oS.K[i, i] - oS.K[j, j]
        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.K[i, i] - oS.labelMat[j] * (
                    oS.alphas[j] - alphaJold) * oS.X[i, j]
        b2 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (
                    oS.alphas[j] - alphaJold) * oS.X[j, j]
        # 步骤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

5.3 在测试中使用核函数

def testRbf(k1=1.3):
    dataArr, labelArr = loadDataSet('D:\\learning\\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('D:\\learning\\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))

jgq1466693

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
机器学习第六章笔记——支持向量机

元组的第一个参数是描述所用核函数类型的一个字符串，其他参数是核函数可能需要的可选参数。在线性核函数的情况下，内积计算在”所有数据集“和数据集中的一行这两个输入之间展开，径向基函数的情况下，在for循环中对于矩阵的每一个元素计算高斯函数的值，当for循环结束后，将计算过程应用到整个向量上去。通过上面观察上面两组数据集，可以看出对于线性可分数据集，可以轻易的画出一条直线将其分开，这条直线又叫分隔超平面，因为给出的例子是二维的，所以此时的超平面就是一条直线，若数据是三维的，那么此时用来分割数据的就是一个平面了。
复制链接

扫一扫