机器学习——支持向量机

定义：

支持向量机：英文全程“ Support Vector Machines ”，简称“ SVM ”，是机器学习中最常用的一种“分类算法”，主要用于构建二分类模型解决二分类问题。

支持向量机是一类按监督学习方式对数据进行二元分类的广义线性分类器，其决策边界是对学习样本求解的最大边距超平面，可以将问题化为一个求解凸二次规划的问题。他可以分为以下两种情况：

• 线性可分：在原空间寻找两类样本的最优分类超平面。
• 线性不可分：加入松弛变量并通过使用非线性映射将低维度输入空间的样本映射到高维度空间使其变为线性可分，这样就可以在该特征空间中寻找最优分类超平面。

首先，先了解以下“支持向量机”是什么意思：

以象棋为例：中国象棋分为黑棋和红棋，并用“ 楚河汉界 ”将其分开。如果用一条直线将不同颜色的棋子进行分类，这显然信手拈来，只需要在楚河汉界的空白附带画一条“中轴线”就能以最佳的方式将它们分开，这样就能保证两边距离最近的棋子保有充分的“间隔”。“ 间隔 ”的产生实际上是依据两侧不同颜色的棋子划分而成的，把这些棋子统称为“样本点”。虽然这些样本点都参与了分类，但对于分类效果影响最大的是处于间隔“边缘”的样本，只要将处于边缘的样本正确分类，那么这两个类别也就分开了，因此我们把处于边缘样本点称为**“支持向量”**。而” 机 “指的是” 一种算法 “。

对于支持向量机而言，它有着重要的三个构件：最大间隔、高维映射、核函数。接下来就一起刨析支持向量机的基本概念吧。

基础概念

最大间隔与分类

线性可分

定义：可以用一个线性函数把两类样本分开，可以是二维平面中的直线，三维空间的平面，多维空间的超平面。

数学定义： 和 是 维欧氏空间中的两个点集。如果存在 维向量 和实数 ，使得所有属于 的点 都有 ，而对于所有属于 的点 则有 ，则我们称 和 线性可分。

以二维图像为例：

我们将数据集分隔开来的直线称为分隔超平面，即。

不同维度的分割超平面：

我们通过观察图形可以发现，超平面将空间分成两部分：

• 点在超平面的“上方” ，记为满足以下式子：

$${\mathbf{w}}^T{\mathbf{x}}_{+}+b>0$$

• 点在超平面的“下方”，记为满足以下式子：
  $${\mathbf{w}}^T{\mathbf{x}}_{-}+b<0$$

• 注意：“上方”和“下方”并不是方位上的超平面上下方，而是以超平面的法向量的指向为准，指向的方向称为“上方”，反之则为“下方”。

• 正负样本对应的标记值为:
  $$y_{+}=+1,y_{-}=-1$$
  所以上述两个式子可以写成：
  $$y_{+}({\mathbf{w}}^T{\mathbf{x}}_{+}+b)>0$$

  $$y_{-}({\mathbf{w}}^T{\mathbf{x}}_{-}+b)>0$$

  故对于线性可分的样本集：

$$D=((x^i,y^i)|i=1,2,...,m),x^i\in R^{1\times n},y^i\in (-1,+1)$$

分类正确的超平面需要满足的条件为:
$$y^i({\mathbf{w}}^T{\mathbf{x}}^{\mathbf{i}}+b)>0$$

$$\widehat{{\gamma }^i}=y^i({\mathbf{w}}^T{\mathbf{x}}^{\mathbf{i}}+b)>0$$

$$\widehat{{\gamma }^i}>0,i=1,2,...m$$

间隔

假设两类数据可以被 （w是超平面的法向量，b是超平面的偏置项）分离，垂直于法向量，移动直到碰到某个训练点，可以得到两个超平面 和 ，两个平面称为支撑超平面，它们分别支撑两类数据。而位于 和 正中间的超平面是分离这两类数据最好的选择。

法向量 有很多中选择，超平面 和 之间的距离称为间隔，这个间隔是 和b的函数，目的就是寻找这样的 和b使得间隔达到最大。

在svm中，将函数间隔亦可以用下面式子表示：
$$\widehat{{\gamma }^i}=y^i({\mathbf{w}}^T{\mathbf{x}}^{\mathbf{i}}+b）$$
由上线性可分中可以了解到，分类正确的超平面需要满足的条件为：样本点到超平面的函数间隔大于零。

硬间隔

定义：决策边界中训练集不能存在分类错误的情况，这一般只存在与理想状态下，现实中很难达到，如下例图。

软间隔

定义：决策边界中训练集存在较小的训练误差，如下例图。

最大间隔

通过上述两个间隔（软间隔和硬间隔）图形中，我们可以了解到无论是软间隔还是硬间隔，我们都可以发现他们的超平面到两个样本的支持向量都存在着“充分”间隔，这个间隔就是最大间隔。为啥要使得两个样本的支持向量保持最大间隔呢？

原因： 如果将数据样本分割的不留余地，即间隔不是最大，就会对随机扰动的噪点特别敏感，这样就很容易破坏掉之前的分类结果，学术上称为“鲁棒性差”。

那怎么从两个分类样本中使得超平面到两个样本的支持向量最大呢？

从上面描述间隔的时候我们就知道，想要找到最大间隔，需要找到其对应的 和b。

1. 计算某一个样本x_i到超平面的距离，我们可以通过几何距离来求解，如下图和公式。

$$d=\frac{|w^T{x}_i+b|}{||w||}=\frac{y_i(w^T{x}_i+b)}{||w||}$$

$$||w||=\sqrt{w_1^2+...+w_n^2}$$

2. 由上述间隔的定义可知，我们要从样本集中找到的是支持向量到超平面的距离，并且得知样本点到超平面的函数间隔大于零，于是可由下式等效：

$$\exists w,bs.t{W}^T{x}_i+b>0if{y}_i=+1$$

$$W^T{x}_i+b<0if{y}_i=-1$$

$$\iff \exists w,bs.t{y}_i(W^T{x}_i+b)>0$$

$$\iff \exists w,b,cs.t{y}_i(W^T{x}_i+b)\ge candc>0$$

将两边的式子同除c，左边即会形成一个新的W矩阵和偏置项b。
$$\exists w,bs.t{y}_i(W^T{x}_i+b)\ge 1i=1,2,..m$$

3. 由1公式，正、负样本到超平面的距离为：

$$d_{+}=\frac{|w^T{x}_i+b|}{||w||}$$

$$d_{-}=\frac{|w^T{x}_i+b|}{||w||}$$

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9KVbi3uJ-1670729760432)(D:\桌面\支持向量机\图片8.png)]

所以可得间隔为：
$$width=\frac{2}{||w||}$$

4. 求解最大化间隔

由3最后得到的width可知最大化间隔即要找到参数w、b使得width最小，即变成求解一个凸二次规划问题：
$$arg\underset{w,b}{\max }\frac{2}{||w||}$$

$$s.tyi(w^T{x}_i\ge 1,i=1,2...,m)$$

由数学知识我们可知：
$$max(\frac{2}{||w||})\iff min(\frac{||w||}{2})$$
为了更好的求导简化计算，我们将其等价于求解下式的最小值
$$arg\underset{w,b}{\min }\frac{1}{2}||w|{|}^2$$
至此，最大间隔的求解完成。

拉格朗日乘子法与对偶问题

拉格朗日乘子法 (Lagrange multipliers)： 一种寻找多元函数在一组约束下的极值的方法。通过引入拉格朗日乘子，可将有 d 个变量与 k 个约束条件的最优化问题转化为具有 d+k 个变量的无约束优化问题求解。

拉格朗日函数的定义：

等式约束：
$$\underset{x}{\max }f(x)$$

$$s.t.g(x)=0$$

拉格朗日函数：
$$\begin{gathered} L(\mathbf{x},\lambda )=f(x)+\lambda g(x) \\ \lambda 是拉格朗日乘子 \end{gathered}$$
于是原本的约束优化问题转换成等价的无约束优化问题：
$$\underset{x,\lambda }{\min }L(x,\lambda )$$
对拉格朗日函数L求偏导：
$${\nabla }_{x}L=\frac{\alpha L}{\alpha x}=\nabla f+\lambda \nabla g=0$$

$${\nabla }_{\lambda }L=\frac{\alpha L}{\alpha y}=g(x)=0$$

不等式约束的KKT条件
$$\underset{x}{\max }f(x)$$

$$s.t.g(x)\le 0$$

可得拉格朗日函数：
$$L(\mathbf{x},\lambda )=f(x)+\lambda g(x)$$
格朗日乘子法的几何意义即在等式g(x)=0或在不等式约束g(x)≤0下最小化目标函数f(x)，阴影部分代表不等式约束表示的可行域。

• g（x）< 0，最优解 x∗在可行域内，此时不等式约束不起作用，只需求解 f(x) 的极值即可。
• g（x）= 0，最优解 x∗ 在可行域边界上，此时不等式约束退化为等式约束。

其约束范围为不等式，因此可等价转化成Karush-Kuhn-Tucker (KKT)条件:
$${\nabla }_{x}L=\nabla f+\lambda \nabla g=0$$

$$\lambda \ge 0$$

$$g(x)\le 0$$

$$\lambda g(x)=0$$

最大间隔问题的拉格朗日乘法：
$$arg\underset{w,b}{\min }\frac{1}{2}||w|{|}^2$$

$$s.t{y}_i(w^T{x}_i+b)\ge 1i=1,2,...m$$

• 第一步：引入拉格朗日乘子 α_i≥0得到拉格朗日函数

$$L(w,b,\alpha )=\frac{1}{2}||w|{|}^2-\sum _{i=1}^m{\alpha }_i(y_i(w^T{x}_i+b)-1)$$

• 第二步：令L(w,b,α) 对w和b的偏导为零

$$w=\sum _{i=1}^m{\alpha }_i{y}_i{x}_i\sum _{i=1}^m{\alpha }_i{y}_i=0$$

• 第三步：w, b回代到第一步，化简

• 第四步：从而得到对偶问题

$$max\sum _{i=1}^m{\alpha }_i-\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i^T{x}_j$$

$$s.t.\sum _{i=1}^m{\alpha }_i{y}_i=0{\alpha }_i\ge 0,i=1,2...m$$

• 第五步：此式为关于的极大值求解，当求出解之后，求出，有

$$f(x)=w^{T}x+b=\sum _{i=1}^m{\alpha }_i{y}_i{x}_i^T{x}_j+b$$

• 由KKT可得：

$${\alpha }_i\ge 0$$

$$y_{i}f(x_{[}i])\ge 1$$

$${\alpha }_i(y_{i}f(x_i)-1)=0$$

对于不在最大边缘边界上的点：
$$y_{i}f(x_i)>1==>{\alpha }_i=0$$
支持向量机解的稀疏性: 训练完成后, 大部分的训练样本都不需保留, 最终模型仅与支持向量有关.

SMO算法与对偶问题

对于大样本数据集的计算量使用拉格朗日乘子的计算量太大了，于是就有了SMO算法。

原理：每次循环中选择两个α进行优化处理。一旦找到一对合适的α，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个α必须要符合 一定的条件，条件之一就是这两个α必须要在间隔边界之外，而其第二个条件则是这两个α还没有进行过区间化处理或者不在边界上。

$$\underset{a}{\max }\sum _{i=1}^m{\alpha }_i-\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i^T{x}_j$$

$$s.t.\sum _{i=1}^m{\alpha }_i{y}_i=0{\alpha }_i\ge 0,i=1,2...m$$

基本思路：不断执行下面两个步骤直到收敛

• 第一步：选取一对需更新的变量α_i和 α_j.

• 第二步：固定α_i和 α_j以外的参数, 求解对偶问题更新α_i和 α_j.

仅考虑α_i和 α_j时, 对偶问题的约束变为
$${\alpha }_i{y}_i+{\alpha }_j{y}_j=-\sum _{k\ne i,j}{\alpha }_k{y}_k,{\alpha }_i\ge 0,{\alpha }_j\ge 0$$

• 偏移项b：通过支持向量来确定

算法流程：每次选择两个α进行更新

为什么要选取两个α进行更新？

如果选择一个，该变量可以通过其它变量和约束条件联合求得下式
$$\sum _P{i=1}^m{\alpha }_i{y}_i=0$$

假设最优解：
$${\alpha }^{\ast }=({\alpha }_1^{\ast },{\alpha }_2^{\ast },...,{\alpha }_n^{\ast })$$

$$w^{\ast }=\sum _{i=1}^m{\alpha }_i^{\ast }{y}_i{x}_i{x}_j$$

可得：
$$b^{\ast }=y_j-\sum _{i=1}^m{\alpha }_i^{\ast }{y}_i{x}_i{x}_j$$
得到超平面：
$$f（x）=w^{\ast }x+b^{\ast }$$

高维映射

高维映射主要是用来解决“你中我，我中有你”的分类问题的，也就是前面所说的“线性不可分问题”，所谓高维映射就是站在更高的维度来解决低维度的问题。

原理：通过增加一个维度的方法，解决“线性不可分的问题”。

点线面可以构成三维立体图，比如棋子是棋盘上的“点"，“间隔”是棋盘上的一条线，棋盘则是一个“面”，而当我们拍盘而起，棋子飞升就会形成一个多维的立体空间，如下图：

经过高维映射后，二维分布的样本点就变成了三维分布，而那张恰好分开棋子的超平面。

核函数

上述高维映射过程是通过核函数（或称映射函数）来实现的，通过这个函数就可以找到一个三维空间，并确定数据点分布，至于能否保证样本点完全分开，这也是由核函数决定的。

定义：通过某非线性变换 φ( x) ，将输入空间映射到高维特征空间。

作用：由低维度空间向高维度空间作一个映射，使原本线性不可分数据变得在高维度上变得可分，并找到这个分割函数。

常见的核函数：

代码实现

SMO算法的代码实现

1. 数据集展示

from numpy import *
import matplotlib.pyplot as plt
 
# 读取数据
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 showData():
    dataMat, labelMat = loadDataSet(r'D:\桌面\SVM\testSet.txt')         # 加载数据集，标签
    dataArr = array(dataMat)                               # 转换成numPy的数组
    n = shape(dataArr)[0]                                  # 获取数据总数
    xcord1 = []; ycord1 = []                               # 存放正样本
    xcord2 = []; ycord2 = []                               # 存放负样本
    for i in range(n):                                     # 依据数据集的标签来对数据进行分类
        if int(labelMat[i]) == 1:                          # 数据的标签为1，表示为正样本
            xcord1.append(dataArr[i, 0]); ycord1.append(dataArr[i, 1])
        else:                                              # 否则，若数据的标签不为1，表示为负样本
            xcord2.append(dataArr[i, 0]); ycord2.append(dataArr[i, 1])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=15, c='blue')             # 绘制正样本
    ax.scatter(xcord2, ycord2, s=15, c='red', marker='s')  # 绘制负样本
    plt.title('DateSet')                                   # 标题
    plt.xlabel('X1'); plt.ylabel('X2')                     # x，y轴的标签
    plt.show()
    
showData()

2. 简化版SMO算法处理小数据集

算法伪代码：

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

代码：

# 随机选择alpha
def selectJrand(i, m):
    j = i                               # 选择一个不等于i的j
    while (j == i):                     # 只要函数值不等于输入值i，函数就会进行随机选择
        j = int(random.uniform(0, m))
    return j
 
# 修剪alpha
def clipAlpha(aj, H, L):                # 用于调整大于H或小于L的alpha值
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

# 简化版SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)              # 数据矩阵dataMatIn转换为numpy的mat存储
    labelMat = mat(classLabels).transpose()  # 数据标签classLabels转换为numpy的mat存储
    b = 0; m, n = shape(dataMatrix)          # 初始化b参数，统计dataMatrix的维度m*n
    alphas = mat(zeros((m, 1)))              # 初始化alpha参数为0
    iter = 0                                 # 初始化迭代次数0
    while (iter < maxIter):                  # matIter表示最多迭代次数，iter变量达到输入值maxIter时，函数结束运行并退出
        alphaPairsChanged = 0                # 变量alphaPairsChanged用于记录alpha是否已经进行优化
        for i in range(m):
            # 步骤1：计算误差Ei
            fXi = float(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)):
                j = selectJrand(i, m)        # 随机选择另一个与alpha_i成对优化的alpha_j
 
                # 步骤1：计算误差Ej
                fXj = float(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("j not moving enough")
                    continue
 
                # 步骤6：更新alpha_i
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])  # 按与alpha_j相同的方法更新alpha_i
 
                # 步骤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, i, alphaPairsChanged))
 
        # 更新迭代次数
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("迭代次数: %d" % iter)
    return b, alphas

# 计算w值
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

# 绘制数据集以及划分直线
def showDataLine(w, b):
    x, y = loadDataSet('testSet.txt')
    xarr = array(x)
    n = shape(x)[0]
    x1 = []; y1 = []
    x2 = []; y2 = []
    for i in arange(n):
        if int(y[i]) == 1:
            x1.append(xarr[i, 0]);
            y1.append(xarr[i, 1])
        else:
            x2.append(xarr[i, 0]);
            y2.append(xarr[i, 1])
 
    plt.scatter(x1, y1, s=30, c='r', marker='s')
    plt.scatter(x2, y2, s=30, c='g')
 
    # 画出 SVM 分类直线
    xx = arange(0, 10, 0.1)
    # 由分类直线 weights[0] * xx + weights[1] * yy1 + b = 0 易得下式
    yy1 = (-w[0] * xx - b) / w[1]
    # 由分类直线 weights[0] * xx + weights[1] * yy2 + b + 1 = 0 易得下式
    yy2 = (-w[0] * xx - b - 1) / w[1]
    # 由分类直线 weights[0] * xx + weights[1] * yy3 + b - 1 = 0 易得下式
    yy3 = (-w[0] * xx - b + 1) / w[1]
    plt.plot(xx, yy1.T)
    plt.plot(xx, yy2.T)
    plt.plot(xx, yy3.T)
 
    # 画出支持向量点
    for i in range(n):
        if alphas[i] > 0.0:
            plt.scatter(xarr[i, 0], xarr[i, 1], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
 
    plt.xlim((-2, 12))
    plt.ylim((-8, 6))
    plt.show()

# 主函数
if __name__ == '__main__':
    dataMat, labelMat = loadDataSet(r'D:\桌面\SVM\testSet.txt')
    b, alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    w = calcWs(alphas, array(dataMat), labelMat)
    showDataLine(w, b)

运行结果：

3. 利用完整Platt SMO算法加速优化

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

# 读取数据
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
def selectJrand(i, m):
    j = i                               # 选择一个不等于i的j
    while (j == i):                     # 只要函数值不等于输入值i，函数就会进行随机选择
        j = int(random.uniform(0, m))
    return j
 
# 修剪alpha
def clipAlpha(aj, H, L):                # 用于调整大于H或小于L的alpha值
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

# 类
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):  # 使用参数初始化结构
        self.X = dataMatIn                                       # 数据矩阵
        self.labelMat = classLabels                              # 数据标签
        self.C = C                                               # 松弛变量
        self.tol = toler                                         # 容错率
        self.m = shape(dataMatIn)[0]                             # 数据矩阵行数m
        self.alphas = mat(zeros((self.m, 1)))                    # 根据矩阵行数初始化alpha参数为0
        self.b = 0                                               # 初始化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 selectJ(i, oS, Ei):
    maxK = -1; maxDeltaE = 0; Ej = 0     # 初始化
    oS.eCache[i] = [1, Ei]               # 选择给出最大增量E的alpha
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:        # 循环使用有效的Ecache值并找到使delta E最大化的值
            if k == i: continue          # 如果k对于i，不计算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
    else:                                # 在这种情况下（第一次），没有任何有效的eCache值
        j = selectJrand(i, oS.m)         # 随机选择alpha_j的索引值
        Ej = calcEk(oS, j)
    return j, Ej

#  计算Ek并更新误差缓存
def updateEk(oS, k):                     # 任何alpha更改后，更新缓存中的新值
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]
    
# 优化的SMO算法
def innerL(i, oS):
    Ei = calcEk(oS, i)                   # 计算误差Ei
    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)):
        # 使用内循环启发方式选择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("j not moving enough"); 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,kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(), C, toler, kTup)# 初始化
    iter = 0                                                                   # 初始化迭代次数为0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):       # 超过最大迭代次数或者遍历整个数据集都alpha也没有更新,则退出循环
        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 = 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):
            entireSet = True
        print("迭代次数: %d" % iter)
    return oS.b, oS.alphas

# 计算w
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

# 绘制数据集以及划分直线
def showData(w, b):
    x, y = loadDataSet('testSet.txt')
    xarr = array(x)
    n = shape(x)[0]
    x1 = []; y1 = []
    x2 = []; y2 = []
    for i in arange(n):
        if int(y[i]) == 1:
            x1.append(xarr[i, 0]);
            y1.append(xarr[i, 1])
        else:
            x2.append(xarr[i, 0]);
            y2.append(xarr[i, 1])
 
    plt.scatter(x1, y1, s=30, c='r', marker='s')
    plt.scatter(x2, y2, s=30, c='g')
 
    # 画出 SVM 分类直线
    xx = arange(0, 10, 0.1)
    # 由分类直线 weights[0] * xx + weights[1] * yy1 + b = 0 易得下式
    yy1 = (-w[0] * xx - b) / w[1]
    # 由分类直线 weights[0] * xx + weights[1] * yy2 + b + 1 = 0 易得下式
    yy2 = (-w[0] * xx - b - 1) / w[1]
    # 由分类直线 weights[0] * xx + weights[1] * yy3 + b - 1 = 0 易得下式
    yy3 = (-w[0] * xx - b + 1) / w[1]
    plt.plot(xx, yy1.T)
    plt.plot(xx, yy2.T)
    plt.plot(xx, yy3.T)
 
    # 画出支持向量点
    for i in range(n):
        if alphas[i] > 0.0:
            plt.scatter(xarr[i, 0], xarr[i, 1], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
 
    plt.xlim((-2, 12))
    plt.ylim((-8, 6))
    plt.show()

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet('testSet.txt')
    b, alphas = smoP(dataMat, labelMat, 0.6, 0.001, 40)
    w = calcWs(alphas, array(dataMat), labelMat)
    showData(w, b)

运行结果：

使用SVM实现手写体数字识别

数据集均来自于——机器学习实战中提供的数据集。

from numpy import *
 
# 随机选择alpha
def selectJrand(i, m):
    j = i                               # 选择一个不等于i的j
    while (j == i):                     # 只要函数值不等于输入值i，函数就会进行随机选择
        j = int(random.uniform(0, m))
    return j
 
# 修剪alpha
def clipAlpha(aj, H, L):                # 用于调整大于H或小于L的alpha值
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj
 
# 类
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):  # 使用参数初始化结构
        self.X = dataMatIn                                       # 数据矩阵
        self.labelMat = classLabels                              # 数据标签
        self.C = C                                               # 松弛变量
        self.tol = toler                                         # 容错率
        self.m = shape(dataMatIn)[0]                             # 数据矩阵行数m
        self.alphas = mat(zeros((self.m, 1)))                    # 根据矩阵行数初始化alpha参数为0
        self.b = 0                                               # 初始化b参数为0
        self.eCache = mat(zeros((self.m, 2)))                    # 第一列是有效标志
        self.K = mat(zeros((self.m,self.m)))                     # 初始化核K
        for i in range(self.m):                                  # 计算所有数据的核K
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
 
# 通过核函数将数据转换更高维的空间
def kernelTrans(X, A, kTup):
    m,n = shape(X)
    K = mat(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 = exp(K/(-1*kTup[1]**2))                      #计算高斯核K
    else: raise NameError('核函数无法识别')
    return K
 
# 计算误差
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
 
# 内循环启发方式
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]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:        # 循环使用有效的Ecache值并找到使delta E最大化的值
            if k == i: continue          # 如果k对于i，不计算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
    else:                                # 在这种情况下（第一次），没有任何有效的eCache值
        j = selectJrand(i, oS.m)         # 随机选择alpha_j的索引值
        Ej = calcEk(oS, j)
    return j, Ej
 
#  计算Ek并更新误差缓存
def updateEk(oS, k):                     # 任何alpha更改后，更新缓存中的新值
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]
 
# 优化的SMO算法
def innerL(i, oS):
    Ei = calcEk(oS, i)                   # 计算误差Ei
    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)):
        # 使用内循环启发方式选择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("j not moving enough"); 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.K[i, j]
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.K[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
 
# 完整的线性SMO算法
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(), C, toler, kTup)# 初始化
    iter = 0                                                                   # 初始化迭代次数为0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):       # 超过最大迭代次数或者遍历整个数据集都alpha也没有更新,则退出循环
        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 = 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):
            entireSet = True
        print("迭代次数: %d" % iter)
    return oS.b, oS.alphas
 
# 图像转换为向量
def img2vector(filename):
    returnVect = zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32 * i + j] = int(lineStr[j])
    return returnVect
 
# 加载图像数据
def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)     # 加载训练集
    m = len(trainingFileList)
    trainingMat = 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
 
# 测试
def testDigits(kTup=('rbf', 10)):
    dataArr, labelArr = loadImages(r'D:\桌面\SVM\digits\trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = mat(dataArr);
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd];
    print("支持向量机是 %d " % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("训练集错误率: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = mat(dataArr);
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("测试错误率: %f" % (float(errorCount) / m))
    
    
if __name__ == '__main__':
    testDigits()

运行结果如下：

总结

支持向量机的优点：

• 有严格的数学理论支持，可解释性强，不依靠统计方法，从而简化了通常的分类和回归问题；
• 能找出对任务至关重要的关键样本（即：支持向量）；
• 采用核技巧之后，可以处理非线性分类/回归任务；
• 最终决策函数只由少数的支持向量所确定，计算的复杂性取决于支持向量的数目，而不是样本空间的维数，这在某种意义上避免了“维数灾难”。

缺点：

• 训练时间长。当采用 SMO 算法时，由于每次都需要挑选一对参数，因此时间复杂度为 O(N2) (N 为训练样本的数量)；

• 当采用核技巧时，如果需要存储核矩阵，则空间复杂度为 O(N2) ；

• 模型预测时，预测时间与支持向量的个数成正比。当支持向量的数量较大时，预测计算复杂度较高。

• 支持向量机目前只适合小批量样本的任务，无法适应百万甚至上亿样本的任务。

代码链接：

链接：https://pan.baidu.com/s/1ZNBE7EZCFSVQwdGqGr1Q2w
提取码：ahsg