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

引入

  有人认为，支持向量机(Support Vector Machines, SVM)是最好的现成的分类器，这里所说的“现成”是指分类器不加修改即可直接使用。同时，这就意味着在数据集上应用基本形式的SVM分类器就可以得到低错误率的结果。SVM能够对训练集之外的数据点作出很好的分类决策。
　　SVM有多种实现，这里只关注其中最流行的一种，即序列最小优化(Sequential Minimal Optimization, SMO)算法，并将使用核函数(kernel)将SVM扩展到更多数据集上。

1 基于最大间隔分隔数据

支持向量机

  优点：泛化错误率低，计算开销不大，结果易解释。
　　缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于二分类问题。
　　适用数据类型：数值型和标称型数据。

  如下图1-1所示，共给定四个数据集：

图1-1　４个示例数据集

 
  请思考这样一个问题：对于A、B、C、D四个数据集中的任意一个，能否只用一条直线将其分开(不必纠结个别点划分错误)？显然Ａ、B数据集时可以的，这组数据便被称为线性可分(linearly separable)，反之则为线性不可分(Linear inseparability)。
 
　　上述将数据集分隔开来的直线称为分隔超平面(separating hyperplane)。以上例子中，所有的数据点都位于二维平面，所有分隔超平面只是一条直线。对于$N$维的情况，则需要$N-1$的某某对象来进行分隔，该对象被称为超平面(hyperplane)，也就是分类的决策边界，作用类似于logistic回归的最佳拟合直线。
 
　　希望构建分类器的方式如下：
　　如果数据点离决策点边界越远，那么最后的预测也就越可信。考虑图1-2(来源)框B到框D的三条直线：
　　

图1-2　一个线性可分的数据集及对应的三种划分情况

 
  B、C、D三个框中直线均能将数据分开，但是哪条最好呢？是否应该最小化数据点到分隔超平面的平均距离来求最佳直线？如果如此，B框和C框中的直线就比D框中的好？这样做是不是有点寻找最佳拟合直线的感激？
　　是的，以上做法确实有点像直线拟合，但这并非最佳方案。我们希望找到离分隔超平面最近的点，以确保它们离分隔面的距离尽可能的远。这里点到分隔面的距离被称为间隔(margin)。我们希望间隔尽可能大，这是因为如果犯错或在有限数据集上训练分类器，分类器应该尽可能健壮。
　　
　　支持向量(support vector)：离分隔超平面最近的点的集合。

2 寻找最大间隔

  如何求解数据集的最佳分隔直线？首先有如下数据集。分隔超平面的形式可写作$W^{T}X+b$。要计算点$A$到分隔超平面的距离，就必须给出点到分隔面的法线或者垂线的长度，该值为$|W^{T}A+b|/||W||$。这里的常数$b$类似于logistic回归中的截距$W_0$。这里的向量$W$和常数$b$一起描述了分隔面或超平面。

图2-1　点A到分隔超平面的距离就是该点到分隔面的法线长度

2.1 分类器求解的优化问题

  理解分类器的工作原理有助于理解基于优化问题的分类器的求解过程。输入数据到类似于$sigmoid$函数的函数，即$f(W^{T}X+b)=f(u)$，当$u<0$时输出-1，反之输出+1。
　　这里的标签为何采用-1和+1，而不是0和1？这是因为-1和+1仅仅相差一个符号，方便数学上的处理。可以通过一个统一的公式来表示间隔或者数据点到分隔超平面的距离，同时不必担心数据到底属于-1类还是+1类。
　　当计算数据点到分隔面的距离并确定分隔面的放置位置时，间隔通过$label\ast (W^{T}X+b)$¹来计算，这就体现了-1和+1的好处。如果数据点处于正方向即+1类，并且离分隔超平面足够远，$W^{T}X+b$以及$label\ast (W^{T}X+b)$将是一个很大的正数，-1类时依然如此。

  现在的目标就是找出分类器定义中的$W$和$b$。为此，必须找到具有最小间隔的数据点，即之前所述的支持向量。一旦找到，就需要对该间隔最大化。记作：
a r g max ⁡ W , b { min ⁡ n ( l a b e l ⋅ ( W T x + b ) ⋅ 1 ∣ ∣ W ∣ ∣ ) } (2-1) arg\max_{W,b} \{ \min_n (label·(W^Tx+b)·\frac{1}{||W||}) \} \tag{2-1} argW,bmax​{nmin​(label⋅(WTx+b)⋅∣∣W∣∣1​)}(2-1)  直接求解2-1式相当困难，所以将其转化为另一种更易求解的形式。首先考虑大括号内的部分：由于对乘积优化是一件很讨厌的事情，因此固定其中一个因子而最大化其他因子。
　　如果令所有的支持向量$label\ast (W^{T}X+b)=1$，那么就可以通过求$\frac{1}{||W||}$的最大值来得到最终解。但是，并非所有数据点的$label\ast (W^{T}X+b)$都等于1，只有那些离分隔超平面最近的点得到的值才是1。离超平面越远，意味着$label\ast (W^{T}X+b)$也越大。

  上述优化问题中，给定了一些约束条件然后求解最优值，因此该问题是一个带约束条件的优化问题。这里的约束条件就是$label\ast (W^{T}X+b)\ge 1.0$.对于此类问题，有一个非常著名的求解方法，即拉格朗日乘子法。通过引入拉格朗日乘子，就可以基于约束条件来表述原来的问题。
　　由于这里的约束条件都是基于数据的的，因此可以将超平面写成数据点的形式。于是最终的优化目标函数可以写成：
$$\underset{\alpha }{\max }\left[\sum _{i=1}^m\alpha -\frac{1}{2}\sum _{i,j=1}^{m}labe{l}^{(i)}\cdot labe{l}^{(j)}\cdot {\alpha }_i\cdot {\alpha }_j<x^{(i)},x^{(j)}>\right]\text{\hspace{1em}(2-2)}$$其中$<x^{(i)},x^{(j)}>$表示$x^{(i)}$和$x^{(j)}$的内积，其约束条件为：
$$\alpha \ge 0，和\sum _{i-1}^m{\alpha }_i\cdot labe{l}^{(i)}=0\text{\hspace{1em}(2-3)}$$  至此，一切都很完美，但是这里有个假设：数据必须100%线性可分。但是，几乎没有这样完美的数据集。这时就可以通过引入松弛变量(slack variable)，来允许有些数据点可以处于分隔面错误一侧。这样优化目标不变，仅仅是约束条件变为：
$$C\ge \alpha \ge 0，和\sum _{i-1}^m{\alpha }_i\cdot labe{l}^{(i)}=0\text{\hspace{1em}(2-4)}$$  这里的C用于控制“最大化间隔”和“保证大部分带你的函数间隔小于1.0”这两个目标的权重。在优化算法的实现代码中，常数C是一个参数，因此可以通过调节C来得到不同的结果。一旦求解出所有$\alpha$，那么分隔超平面就可以通过这些$\alpha$来表示。这一结论十分直接，SVM的主要工作就是求解$\alpha$。

2.2 SVM的一般框架

  (1)收集数据：任意方法；
　　(2)准备数据：数值型数据；
　　(3)分析数据：有助于可视化分隔超平面；
　　(4)训练算法：主要实现两个参数的调优；
　　(5)测试算法：十分简单的计算过程即可实现；
　　(6)使用算法：几乎所有的分类问题都可以使用SVM，值得一提的是，SVM本身是一个二分类分类器，对于多类问题需要适当修改。

3 SMO高效优化算法

  优化的目标有二：
　　(1)最小化目标函数；
　　(2)优化过程中必须遵循的约束条件。
　　所有需要围绕优化做的事就是训练分类器，一旦得到$\alpha$的最优值，就得到了分隔超平面并用之于分类。

3.1 Platt的SMO算法

  1996年，John Platt发布了一个称为SMO²的强大算法，用于训练SVM。SMO表示序列最小优化(Sequential Minimal Optimization)。Platt的SMO算法是将大优化问题分解为小优化问题来求解。这些小优化问题往往容易求解，并且它们进行顺序求解的结果将于=与它们作为整体来求解的结果完全一致，并能大幅降低求解时间。
 
　　SMO算法的工作原理：
　　每次循环中选择两个$\alpha$进行优化处理。一旦找到了一对合适的$\alpha$，那么就增加其中一个的同时减小另一个。这里的“合适”就是指两个$\alpha$必须要符合一定的条件：1）两个$\alpha$必须要在间隔边界之外；2）两个$\alpha$还没有进行过区间化处理或者不在边界上。

3.2 简化版的SMO算法

  简化版的SMO算法代码量少，方便了解算法原理，不过运行速度慢，最好用于小规模数据集。
　　Platt SMO算法中的外循环确定要优化的最佳$\alpha$对。而简化版会跳过这一步骤，首先在数据集上遍历每一个$\alpha$，然后在剩下的$\alpha$集合中随机选择另一个$\alpha$，从而构建$\alpha$对。这里有一点相当重要，就是需要同时改变两个$\alpha$。之所以这么做是因为必须满足以下约束条件：
$$\sum {\alpha }_i\cdot labe{l}^{(i)}=0\text{\hspace{1em}(3-1)}$$  由于改变一个$\alpha$会导致约束条件失效，因此总是同时改变两个$\alpha$。
　　为此，将构建两个辅助函数，其一用于在某个范围内随机选取一个整数，其二用于在数值过大时进行调整。创建svm_MLiA.py文件并添加以下代码：

**程序清单3-1：**SMO算法的辅助函数

import random

def load_data_set(file_name):
    with open(file_name) as fd:
        fd_data = fd.readlines()

    data_set = []; label_set = []
    for data in fd_data:
        data = data.strip().split('\t')
        data = [float(value) for value in data]
        data_set.append(data[: -1])
        label_set.append(data[-1])
    return data_set, label_set

def select_j_rand(i, m):    #选择一个不等于i的j
    j = i
    while j == i:
        j = int(random.uniform(0, m))
    return j

def clip_alpha(aj, H, L):    #数值调整
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

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

  于svm_MLiA.py文件添加以下代码：
　　　　
**程序清单3-2：**简化版SMO算法

def smo_simple(data_set, label_set, C, toler, max_iter):    #输入参数：数据集、标签集、常数C、容错率、最大循环次数
    data_mat = mat(data_set)
    label_mat = mat(label_set).T
    b = 0
    m, n = shape(data_mat)
    alphas = mat(zeros((m, 1)))    #初始化alphas
    iter = 0    #初始化循环变量
    while iter < max_iter:    #只有在数据集上遍历amx_iter次，且不再发生任何alpha修改之后才会退出循环
        alpha_pairs_changed = 0    #记录alpha是否已经优化
        for i in range(m):
            f_Xi = float(multiply(alphas, label_mat).T * (data_mat * data_mat[i,:].T)) + b     #预测类别
            Ei = f_Xi - float(label_mat[i])    #误差
            if ((label_mat[i] * Ei < -toler) and (alphas[i] < C)) or ((label_mat[i] * Ei > toler) and (alphas[i] > 0)):   #如果误差很大，则优化
                j = select_j_rand(i, m)    #随机选择第二个alpha值
                f_Xj = float(multiply(alphas, label_mat).T * (data_mat * data_mat[j,:].T)) + b    #预测类别
                Ej = f_Xj - float(label_mat[j])    #计算误差
                alpha_iold = alphas[i].copy()    #复制
                alpha_jold = alphas[j].copy()    #复制
                if label_mat[i] != label_mat[j]:    #L和H均用于调整alpha[j]，使其＞0小于C
                    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;
                eta = 2.0 * data_mat[i,:] * data_mat[j,:].T - data_mat[i,:] * data_mat[i,:].T - \
                      data_mat[j,:] * data_mat[j,:].T    #eta为alpha[j]的最优修改量
                if eta >= 0:    #由于eta=0不常发生且计算复杂，故若eta=0时，则直接忽略
                    print("eta >= 0")
                    continue
                alphas[j] -= label_mat[j] * (Ei - Ej) / eta    #新计算alpha[j]
                alphas[j] = clip_alpha(alphas[j], H, L)    #调整
                if abs(alphas[j] - alpha_jold) < e-5:    #检查alpha[j]是否出现轻微变化，若是则退出循环
                    print("j not moving enough")
                    continue
                alphas[i] += label_mat[j] * label_mat[i] * (alpha_jold - alphas[j])    #alpha[i]与alpha[j]变化同样大小，但是方向相反
                b1 = b - Ei - label_mat[i] * (alphas[i] - alpha_iold) * data_mat[i,:] * data_mat[i,:].T - \
                    label_mat[j] * (alphas[j] - alpha_jold) * data_mat[i,:] * data_mat[j,:].T
                b2 = b - Ej - label_mat[i] * (alphas[i] - alpha_iold) * data_mat[i,:] * data_mat[j,:].T - \
                    label_mat[j] * (alphas[j] - alpha_jold) * data_mat[j,:] * data_mat[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
                alpha_pairs_changed += 1    #若未执行任意一条continue，则alpha_pairs_changed+1
                print("iter: %d i: %d, pairs changed %d" % (iter, i ,alpha_pairs_changed))
        if alpha_pairs_changed == 0:    #检测alpha_pairs_changed是否变化
            iter += 1
        else:
            iter = 0
        print("iteration number:", iter)
    return b, alphas

def test1():
    data_set, label_set = load_data_set('testSet.txt')
    b, alphas = smo_simple(data_set, label_set, 0.6, 0.001, 40)
    print("b:", b)
    print("alphas > 0 is:", alphas[alphas > 0])
    print("The support vectoc is:")
    for i in range(len(data_set)):
        if alphas[i] > 0.0:
            print(data_set[i], label_set[i])

if __name__ == '__main__':
    test1()

  运行结果：

L == H
L == H
L == H
iter: 0 i: 3, pairs changed 1
L == H
...
iteration number: 0
iter: 0 i: 17, pairs changed 1
iter: 0 i: 29, pairs changed 2
iteration number: 0
iteration number: 1
...
iteration number: 38
iteration number: 39
iteration number: 40
b: [[-3.8387274]]
alphas > 0 is: [[0.12742924 0.24143342 0.36886266]]
The support vectoc is:
[4.658191, 3.507396] -1.0
[3.457096, -0.082216] -1.0
[6.080573, 0.418886] 1.0

  在数据集上绘制出支持向量：

图3-1　简化版SMO运行结果

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

  简化版SMO算法一个很大的问题便在于运行速度上，如果仅仅用于小数据集，是没有问题的，但是在更大的数据集上运行速度就会变慢。
　　完整版SMO算法实现alpha更改、代数运算的优化环节与简化版SMO一模一样，唯一不同的就是选择alpha的方式，应用了一些能够快速提速的启发方法。
　　
　　Platt SMO算法是通过一个外循环来选择第一个alpha值，并且其选择过程会在两种方式之间交替：
　　1）在所有数据集上进行单遍扫描；
　　2）在非边界alpha中实现单遍扫描。
　　所谓非边界alpha是指：不等于边界0或C的alpha值。对整个数据集的扫描相等容易，而实现非边界alpha的扫描时，首先需要建立这些alpha的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。
　　第一个alpha值选择之后，算法会通过一个内循环选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中，会在选择j之后计算错误率Ej。但在这里，会建立一个全局变量用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。
　　改进算法之前，有必要对之前所有代码进行清理。以下程序清单中包含1个用于清理代码的数据结构和3个用于对E进行缓存的辅助函数。于svm_MLiA.py文件添加以下代码：

**程序清单4-1：**完整版Platt SMO算法辅助函数

class optStruct:    #建立数据结果保存所有值
    def __init__(self, data_set, label_set, C, toler):
        self.data_mat = mat(data_set)
        self.label_mat = mat(label_set)
        self.C = C
        self.toler = toler
        self.m = shape(self.data_mat)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.e_cache = mat(zeros(self.m, 2))    #e_cache第一列给出e_cache是否有效的标志位；第二列给出实际E值

"""计算E值并返回"""
def calc_Ek(o_S, k):    #传入参数：数据结构、索引
    f_Xk = float(multiply(o_S.alphas, o_S.label_mat).T * (o_S.data_mat * o_S.data_mat[k,:].T)) + o_S.b
    Ek = f_Xk - float(o_S.label_mat[k])
    return Ek

"""选择第二个alpha值"""
def select_j(i, o_S, Ei):    #传入参数：索引i、数据结构、误差Ei
    max_K = -1; max_delta_E = 0; Ej = 0
    o_S.e_cache[i] = [1, Ei]    #将输入值Ei设置为“有效的”；有效的是指已经计算完毕
    valid_e_cache_list = nonzero(o_S.e_cache[:,0].A)[0]    #nonzero()返回e_cache[:,0].A中非零值的索引；.A表示将矩阵转化为数组
    if len(valid_e_cache_list) > 1:
        for k in valid_e_cache_list:
            if k == i:
                continue
            Ek = calc_Ek(o_S, k)
            delta_E = abs(Ei - Ek)
            if delta_E > max_delta_E:
                max_K = k; max_delta_E = delta_E; Ej = Ek
        return max_K, Ej
    else:
        j = select_j_rand(i, o_S.m)
        Ej = calc_Ek(o_S, j)
    return j, Ej

"""更新误差值"""
def update_Ek(o_S, k):    #传入参数：数据结构、索引k
    Ek = calc_Ek(o_S, k)
    o_S.e_cache[k] = [1, Ek]

  以下代码配合辅助函数，便可用于寻找决策边界。于svm_MLiA.py文件添加以下代码：

**程序清单4-2：**完整版Platt SMO中的优化例程

"""与简化版SMO几乎一样，除了使用数据结构作为输入、使用select_j更新j，且在alpha改变时更新e_cache"""
def inner_L(i, o_S):    #输入参数：索引、数据结构--与简化版SMO不同之处1
    Ei = calc_Ek(o_S, i)
    if ((o_S.label_mat[i] * Ei< -o_S.toler) and (o_S.alphas[i] < o_S.C)) or \
        ((o_S.label_mat[i] * Ei > o_S.toler) and (o_S.alphas[i] > 0)):
        j, Ej = select_j(i, o_S, Ei)    #与简化版SMO不同之处2：使用select_j而不是select_j_rand更新j
        alpha_iold = o_S.alphas[i].copy()
        alpha_jold = o_S.alphas[j].copy()
        if (o_S.label_mat[i] != o_S.label_mat[j]):
            L = max(0, o_S.alphas[j] - o_S.alphas[i])
            H = min(o_S.C, o_S.C + o_S.alphas[j] - o_S.alphas[i])
        else:
            L = max(0, o_S.alphas[j] + o_S.alphas[i] - o_S.C)
            H = min(o_S.C, o_S.alphas[j] + o_S.alphas[i])
        if L == H:
            print("L == H")
            return 0
        eta = 2.0 * o_S.data_mat[i,:] * o_S.data_mat[j,:].T - o_S.data_mat[i,:] * o_S.data_mat[i,:].T - \
              o_S.data_mat[j, :] * o_S.data_mat[j, :].T
        if (eta >= 0):
            print("eta >= 0")
            return 0
        o_S.alphas[j] -= o_S.label_mat[j] * (Ei - Ej) / eta
        o_S.alphas[j] = clip_alpha(o_S.alphas[j], H, L)
        update_Ek(o_S, j)    #与简化版SMO不同之处3
        if (abs(o_S.alphas[j] - alpha_jold) < e-5):
            print("j not moving enough")
            return 0
        o_S.alphas[i] += o_S.label_mat[j] * o_S.label_mat[i] * (alpha_jold - o_S.alphas)
        update_Ek(o_S, i)
        b1 = o_S.b - Ei - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * \
            o_S.data_mat[i,:] * o_S.data_mat[i,:].T - o_S.label_mat[j] * \
            (o_S.alphas[j] - alpha_jold) * o_S.data_mat[i,:]*o_S.data_mat[j,:].T
        b2 = o_S.b - Ej - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * \
            o_S.data_mat[i, :] * o_S.data_mat[j, :].T - o_S.label_mat[j] * \
            (o_S.alphas[j] - alpha_jold) * o_S.data_mat[j,] * o_S.data_mat[j, :].T
        if (0 < o_S.alphas[i]) and (o_S.C > o_S.alphas[i]):
            o_S.b = b1
        elif (0 < o_S.alphas[j]) and (o_S.C > o_S.alphas[j]):
            o_S.b = b2
        else:
            o_S.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0

  以下代码将上述过程打包，即选择第一个alpha的外循环。于svm_MLiA.py文件添加以下代码：

**程序清单4-3：**完整版Platt SMO的外循环

"""完整版SMO算法"""
def smop(data_set, label_set, C, toler, max_iter, k_tup=('lin', 0):    #输入参数：数据集、标签集、常数C、容错率、最大循环次数、核函数选择(暂时不用)
    o_S = optStruct(data_set, label_set, C, toler)    #初始化数据结构
    iter = 0
    entire_set = True    #是否遍历完整个数据集
    alpha_pairs_changed = 0
    while (iter < max_iter) and ((alpha_pairs_changed > 0) or (entire_set)):
        alpha_pairs_changed = 0
        if entire_set:
            for i in range(o_S.m):
                alpha_pairs_changed += inner_L(i, o_S)
                print("full set, iter: %d i: %d, pairs changed %d" % (iter, i, alpha_pairs_changed))
        else:
            non_bound_is = nonzero((o_S.alphas.A > 0) * (o_S.alphas.A < C))[0]
            for i in non_bound_is:
                alpha_pairs_changed += inner_L(i, o_S)
                print("non-bound, iter: %d i: %d, pairs changed %d" % (iter, i, alpha_pairs_changed))
        iter += 1    #此处的迭代定义为一次循环过程
        if entire_set:
            entire_set = False
        elif (alpha_pairs_changed == 0):
            entire_set = True
        print("iteration number: %d" % iter)
    return o_S.b, o_S.alphas

"""计算W"""
def calc_ws(alphas, data_set, label_set):
    data_mat = mat(data_set)
    label_mat = mat(label_set).T
    m , n =shape(data_mat)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * label_mat[i], data_mat[i,:].T)
    return w

def test2():
    data_set, label_set = load_data_set('testSet.txt')
    b, alphas = smop(data_set, label_set, 0.6, 0.001, 40)
    print("b:",b)
    print("alphas:",alphas[alphas > 0])
    ws= calc_ws(alphas, data_set, label_set)
    print("ws:",ws.T)
    print("The support vectoc is:")
    for i in range(len(data_set)):
        if alphas[i] > 0.0:
            print(data_set[i], label_set[i])

if __name__ == '__main__':
    test2()
    #test1()

  运行结果：

...
full set, iter: 2 i: 98, pairs changed 0
full set, iter: 2 i: 99, pairs changed 0
iteration number: 3
b: [[-2.89901748]]
alphas: [[0.06961952 0.0169055  0.0169055  0.0272699  0.04522972 0.0272699
  0.0243898  0.06140181 0.06140181]]
ws: [[ 0.65307162 -0.17196128]]
The support vectoc is:
[3.542485, 1.977398] -1.0
[2.114999, -0.004466] -1.0
[8.127113, 1.274372] 1.0
[4.658191, 3.507396] -1.0
[8.197181, 1.545132] 1.0
[7.40786, -0.121961] 1.0
[6.960661, -0.245353] 1.0
[6.080573, 0.418886] 1.0
[3.107511, 0.758367] -1.0

  绘制如下：

图4-1　完整版SMO算法得到的支持向量

  calc_ws()函数中，最重要的部分就是for循环部分，虽然在循环中实现的仅仅是多个数的乘积。但是在alphas中，大多数的值都是0，而非零alpha所对应的便是支持向量，最终起作用的也是这些支持向量。对于这些0值，对于w的计算毫无帮助，所以数据集中相对于的数据点也就会被很容易地舍弃。
　　有了b和w之后，便可以进行数据集标签的预测。于svm_MLiA.py文件添加以下代码：：
　　
**程序清单4-4：**测试函数

def test3():
    data_set, label_set = load_data_set('testSet.txt')
    b, alphas = smop(data_set, label_set, 0.6, 0.001, 40)
    ws = calc_ws(alphas, data_set, label_set)
    error = 0.0
    data_mat = mat(data_set)
    for i in range(len(data_mat)):
        predicted_label = sign((data_mat[i] * mat(ws) + b))
        if predicted_label != label_set[i]:
            error += 1
    print("The error rate is:", (error / len(data_mat)))

if __name__ == '__main__':
    test3()
    #test2()
    #test1()

  运行结果：

...
The error rate is: 0.0

  目前测试结果好一个很重要的原因在于数据集本身线性可分。但是倘若两类数据点分别分布在一个圆的内部和外部，那会是怎样的一种情况呢？

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

  考虑最开始时图1-1C，这是被我们称之为线性不可分的情况。显然，该数据中存在某种可以识别的模式。其中一个问题就是，能否像线性情况一样，利用强大的工具来捕捉数据中的这种模式？
　　一种工具被称作核函数(kernelk)，用于将数据转换为易于分类器理解的形式。

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

  图1-1C中，数据处于一个圆中，我们自身很容易意识到这一点。然后，对于分类器而言，只能识别最终的结果是大于0还是小于0。如果将数据从一个特征空间转换到另一个特征空间，在新空间下，就可以利用已有工具对数据进行处理。数学家喜欢将这个过程称之为从一个特征空间到另一个特征空间的映射。通常情况下，为低维特征空间映射到高维特征空间。
　　这个映射的过程是通过核函数实现的。可以将其理解为一个包装器(wrapper)或者接口(interface)，它能把数据从某个很难处理的形式转换为另一种较容易处理的形式。简单起见，可以将核函数想象成一种距离计算方法。
　　
　　SVM优化中一个特别好的地方就是：
　　所有的运算都可以写成内积(inner product，也称点积)的形式。向量的内积是指两个向量相乘，之后得到单个标量或者数值。我们可以将内积替换成核函数，该方式被称为核技巧(kernel trick)或者核“变电”(kernel substation)。
　　众多机器学习算法都可以运用核函数，其中一种较流行的，便是径向基核函数。

5.2 径向基核函数

  径向基核函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以是<0,0>向量或者其他向量开始计算的距离。径向基核函数的高斯版本如下：
$$k(x,y)=exp(\frac{-||x-y|{|}^2}{2{\sigma }^2})\text{\hspace{1em}(6-1)}$$其中，$\sigma$是用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数。
　　上述高斯核函数将数据从其特征空间映射到更高维的空间，准确来说是映射到无穷维的空间(目前无需担心何谓无穷维)。
　　对于图1-1C中的数据，数据点基本上都在一个圆内，示例数据集(训练集testSetRBF.txt)如下：

-0.214824	0.662756	-1.000000
-0.061569	-0.091875	1.000000
0.406933	0.648055	-1.000000
0.223650	0.130142	1.000000
0.231317	0.766906	-1.000000
-0.748800	-0.531637	-1.000000
-0.557789	0.375797	-1.000000
0.207123	-0.019463	1.000000
0.286462	0.719470	-1.000000
0.195300	-0.179039	1.000000
-0.152696	-0.153030	1.000000
0.384471	0.653336	-1.000000
-0.117280	-0.153217	1.000000
-0.238076	0.000583	1.000000
-0.413576	0.145681	1.000000
0.490767	-0.680029	-1.000000
0.199894	-0.199381	1.000000
-0.356048	0.537960	-1.000000
-0.392868	-0.125261	1.000000
0.353588	-0.070617	1.000000
0.020984	0.925720	-1.000000
-0.475167	-0.346247	-1.000000
0.074952	0.042783	1.000000
0.394164	-0.058217	1.000000
0.663418	0.436525	-1.000000
0.402158	0.577744	-1.000000
-0.449349	-0.038074	1.000000
0.619080	-0.088188	-1.000000
0.268066	-0.071621	1.000000
-0.015165	0.359326	1.000000
0.539368	-0.374972	-1.000000
-0.319153	0.629673	-1.000000
0.694424	0.641180	-1.000000
0.079522	0.193198	1.000000
0.253289	-0.285861	1.000000
-0.035558	-0.010086	1.000000
-0.403483	0.474466	-1.000000
-0.034312	0.995685	-1.000000
-0.590657	0.438051	-1.000000
-0.098871	-0.023953	1.000000
-0.250001	0.141621	1.000000
-0.012998	0.525985	-1.000000
0.153738	0.491531	-1.000000
0.388215	-0.656567	-1.000000
0.049008	0.013499	1.000000
0.068286	0.392741	1.000000
0.747800	-0.066630	-1.000000
0.004621	-0.042932	1.000000
-0.701600	0.190983	-1.000000
0.055413	-0.024380	1.000000
0.035398	-0.333682	1.000000
0.211795	0.024689	1.000000
-0.045677	0.172907	1.000000
0.595222	0.209570	-1.000000
0.229465	0.250409	1.000000
-0.089293	0.068198	1.000000
0.384300	-0.176570	1.000000
0.834912	-0.110321	-1.000000
-0.307768	0.503038	-1.000000
-0.777063	-0.348066	-1.000000
0.017390	0.152441	1.000000
-0.293382	-0.139778	1.000000
-0.203272	0.286855	1.000000
0.957812	-0.152444	-1.000000
0.004609	-0.070617	1.000000
-0.755431	0.096711	-1.000000
-0.526487	0.547282	-1.000000
-0.246873	0.833713	-1.000000
0.185639	-0.066162	1.000000
0.851934	0.456603	-1.000000
-0.827912	0.117122	-1.000000
0.233512	-0.106274	1.000000
0.583671	-0.709033	-1.000000
-0.487023	0.625140	-1.000000
-0.448939	0.176725	1.000000
0.155907	-0.166371	1.000000
0.334204	0.381237	-1.000000
0.081536	-0.106212	1.000000
0.227222	0.527437	-1.000000
0.759290	0.330720	-1.000000
0.204177	-0.023516	1.000000
0.577939	0.403784	-1.000000
-0.568534	0.442948	-1.000000
-0.011520	0.021165	1.000000
0.875720	0.422476	-1.000000
0.297885	-0.632874	-1.000000
-0.015821	0.031226	1.000000
0.541359	-0.205969	-1.000000
-0.689946	-0.508674	-1.000000
-0.343049	0.841653	-1.000000
0.523902	-0.436156	-1.000000
0.249281	-0.711840	-1.000000
0.193449	0.574598	-1.000000
-0.257542	-0.753885	-1.000000
-0.021605	0.158080	1.000000
0.601559	-0.727041	-1.000000
-0.791603	0.095651	-1.000000
-0.908298	-0.053376	-1.000000
0.122020	0.850966	-1.000000
-0.725568	-0.292022	-1.000000

  绘制如下：

图6-1　训练数据集

测试集testSetRBF2.txt如下：

0.676771	-0.486687	-1.000000
0.008473	0.186070	1.000000
-0.727789	0.594062	-1.000000
0.112367	0.287852	1.000000
0.383633	-0.038068	1.000000
-0.927138	-0.032633	-1.000000
-0.842803	-0.423115	-1.000000
-0.003677	-0.367338	1.000000
0.443211	-0.698469	-1.000000
-0.473835	0.005233	1.000000
0.616741	0.590841	-1.000000
0.557463	-0.373461	-1.000000
-0.498535	-0.223231	-1.000000
-0.246744	0.276413	1.000000
-0.761980	-0.244188	-1.000000
0.641594	-0.479861	-1.000000
-0.659140	0.529830	-1.000000
-0.054873	-0.238900	1.000000
-0.089644	-0.244683	1.000000
-0.431576	-0.481538	-1.000000
-0.099535	0.728679	-1.000000
-0.188428	0.156443	1.000000
0.267051	0.318101	1.000000
0.222114	-0.528887	-1.000000
0.030369	0.113317	1.000000
0.392321	0.026089	1.000000
0.298871	-0.915427	-1.000000
-0.034581	-0.133887	1.000000
0.405956	0.206980	1.000000
0.144902	-0.605762	-1.000000
0.274362	-0.401338	1.000000
0.397998	-0.780144	-1.000000
0.037863	0.155137	1.000000
-0.010363	-0.004170	1.000000
0.506519	0.486619	-1.000000
0.000082	-0.020625	1.000000
0.057761	-0.155140	1.000000
0.027748	-0.553763	-1.000000
-0.413363	-0.746830	-1.000000
0.081500	-0.014264	1.000000
0.047137	-0.491271	1.000000
-0.267459	0.024770	1.000000
-0.148288	-0.532471	-1.000000
-0.225559	-0.201622	1.000000
0.772360	-0.518986	-1.000000
-0.440670	0.688739	-1.000000
0.329064	-0.095349	1.000000
0.970170	-0.010671	-1.000000
-0.689447	-0.318722	-1.000000
-0.465493	-0.227468	-1.000000
-0.049370	0.405711	1.000000
-0.166117	0.274807	1.000000
0.054483	0.012643	1.000000
0.021389	0.076125	1.000000
-0.104404	-0.914042	-1.000000
0.294487	0.440886	-1.000000
0.107915	-0.493703	-1.000000
0.076311	0.438860	1.000000
0.370593	-0.728737	-1.000000
0.409890	0.306851	-1.000000
0.285445	0.474399	-1.000000
-0.870134	-0.161685	-1.000000
-0.654144	-0.675129	-1.000000
0.285278	-0.767310	-1.000000
0.049548	-0.000907	1.000000
0.030014	-0.093265	1.000000
-0.128859	0.278865	1.000000
0.307463	0.085667	1.000000
0.023440	0.298638	1.000000
0.053920	0.235344	1.000000
0.059675	0.533339	-1.000000
0.817125	0.016536	-1.000000
-0.108771	0.477254	1.000000
-0.118106	0.017284	1.000000
0.288339	0.195457	1.000000
0.567309	-0.200203	-1.000000
-0.202446	0.409387	1.000000
-0.330769	-0.240797	1.000000
-0.422377	0.480683	-1.000000
-0.295269	0.326017	1.000000
0.261132	0.046478	1.000000
-0.492244	-0.319998	-1.000000
-0.384419	0.099170	1.000000
0.101882	-0.781145	-1.000000
0.234592	-0.383446	1.000000
-0.020478	-0.901833	-1.000000
0.328449	0.186633	1.000000
-0.150059	-0.409158	1.000000
-0.155876	-0.843413	-1.000000
-0.098134	-0.136786	1.000000
0.110575	-0.197205	1.000000
0.219021	0.054347	1.000000
0.030152	0.251682	1.000000
0.033447	-0.122824	1.000000
-0.686225	-0.020779	-1.000000
-0.911211	-0.262011	-1.000000
0.572557	0.377526	-1.000000
-0.073647	-0.519163	-1.000000
-0.281830	-0.797236	-1.000000
-0.555263	0.126232	-1.000000

  绘制如下：

图6-2　测试数据集

.  对于这个例子，可以直接检查原始数据，并意识到只需要度量数据点到圆心的距离即可。但是其它形式的数据集呢？高斯核函数便能解决这个问题。于svm_MLiA.py文件添加以下代码：
　　
**程序清单5-1：**核转换函数

def kernel_trans(data_set, data_set_line, k_tup):	 #输入参数：特指支持向量、数据集每行、核函数选择元组
    data_mat = mat(data_set)
    m, n = shape(data_mat)
    K = mat(zeros((m, 1)))
    """此处只给出两种核函数的类型，可以继续添加"""
    if k_tup[0] == 'lin':    #线性核函数：内积计算在“所有数据集”和“数据集中的一行中展开”
        K = data_mat * data_set_line.T
    elif k_tup[0] == 'rbf':    #径向基核函数
        for j in range(m):
            delta_row = data_mat[j,:] - data_set_line
            K[j] = delta_row * delta_row.T
        K = exp(K / (-1 * k_tup[1]**2))    #numpy中/代表对元素展开即元素间的除法，而不像在MATLAB中一样计算矩阵的逆
    else:    #遇到无法识别的元组则退出
        raise NameError('Houston We Have a Problem -- That Kernel is not recogized')
    return K

.  为了使用核函数，还需对代码中的其他部分进行修改。于svm_MLiA.py文件修改以下部分代码(注意修改前后)：
　　
**程序清单5-2：**待修改函数及结构体

"""新增k_tup，为包含核函数信息的元组"""
class optStruct:    #建立数据结果保存所有值
    def __init__(self, data_set, label_set, C, toler, k_tup):    
        self.data_mat = mat(data_set)
        self.label_mat = mat(label_set).transpose()
        self.C = C
        self.toler = toler
        self.m = shape(self.data_mat)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.e_cache = mat(zeros((self.m, 2)))    #e_cache第一列给出e_cache是否有效的标志位；第二列给出实际E值
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:,i] = kernel_trans(self.data_mat, self.data_mat[i,:], k_tup)

def inner_L(i, o_S):    #输入参数：索引、数据结构--与简化版SMO不同之处1
   ...
        """修改1"""
        eta = 2.0 * o_S.K[i, j] - o_S.K[i, i] - o_S.K[j, j]
        #eta = 2.0 * o_S.data_mat[i, :] * o_S.data_mat[j, :].T - o_S.data_mat[i, :] * o_S.data_mat[i, :].T - o_S.data_mat[j, :] * o_S.data_mat[j, :].T
        """修改2"""
        b1 = o_S.b - Ei - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.K[i, i] - \
            o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.K[i, j]
        b2 = o_S.b - Ej - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.K[i, i] - \
        o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.K[j, j]
        # b1 = o_S.b - Ei - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.data_mat[i, :] * o_S.data_mat[i, :].T -  \
        #      o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.data_mat[i, :] * o_S.data_mat[j, :].T
        # b2 = o_S.b - Ej - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.data_mat[i, :] * o_S.data_mat[j, :].T - \
        #      o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.data_mat[j, :] * o_S.data_mat[j, :].T

def calc_Ek(o_S, k):    #传入参数：数据结构、索引
    """修改1"""
    f_Xk = float(multiply(o_S.alphas, o_S.label_mat).T * o_S.K[:,k] + o_S.b) 
    #f_Xk = float(multiply(o_S.alphas, o_S.label_mat).T * (o_S.data_mat * o_S.data_mat[k,:].T)) + o_S.b
    ...

"""修改1"""
def smop(data_set, label_set, C, toler, max_iter, k_tup):
    o_S = optStruct(data_set, label_set, C, toler, k_tup) 

5.3 在测试中使用核函数

.  接下来构建的分类器将使用径向基核函数，该分类器只有一个输入参数，即之前所述的$\sigma$。于svm_MLiA.py文件添加以下代码：

**程序清单5-3：**利用核函数进行分类的径向基测试函数

def test_rbf(k1=1.3):
    data_set, label_set = load_data_set('testSetRBF.txt')
    b, alphas = smop(data_set, label_set, 3, 0.0001, 10000, ('rbf', k1))
    data_mat = mat(data_set); label_mat = mat(label_set).T
    sv_ind = nonzero(alphas.A > 0)[0]
    data_sv = data_mat[sv_ind]; label_sv = label_mat[sv_ind]
    print("There are %d Support Vectors" % shape(data_sv)[0])
    m , n = shape(data_mat)
    error_count = 0
    for i in range(m):
        kernel_eval = kernel_trans(data_sv, data_mat[i,:], ('rbf', k1))
        predict = kernel_eval.T * multiply(label_sv, alphas[sv_ind]) + b
        if sign(predict) != sign(label_set[i]):
            error_count += 1
    print("The training error rate is: %f" % (float(error_count) / m))
    error_count = 0
    data_set, label_set = load_data_set('testSetRBF2.txt')
    data_mat = mat(data_set)
    m, n = shape(data_mat)
    for i in range(m):
        kernel_eval = kernel_trans(data_sv, data_mat[i,:], ('rbf', k1))
        predict = kernel_eval.T * multiply(label_sv, alphas[sv_ind]) + b
        if sign(predict) != sign(label_set[i]):
            error_count += 1
    print("The test error rate is: %f" % (float(error_count) / m))

if __name__ == '__main__':
    test_rbf(0.5)
    #test3()
    #test2()
    #test1()

.  运行结果：

...
There are 23 Support Vectors
The training error rate is: 0.010000
The test error rate is: 0.040000

.  通过设置不同的k1值，也会出现不同的错误率，当时支持向量的数目也会相应增减。如果支持向量数目过少，则会得到很差的决策边界，但是如果过多，也就是相当于每次利用所有数据来进行分类，这种方式即k近邻。
　　　至此的分类都是二分类，如果想了解SMO多分类，建议阅读C.E.Huset等人发表的“A Comparison of Methods for Multiclass Support Vector Machines”。

附录：完整代码

import random
from numpy import *

class optStruct:    #建立数据结果保存所有值
    def __init__(self, data_set, label_set, C, toler, k_tup):    #新增k_tup，为包含核函数信息的元组
        self.data_mat = mat(data_set)
        self.label_mat = mat(label_set).transpose()
        self.C = C
        self.toler = toler
        self.m = shape(self.data_mat)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.e_cache = mat(zeros((self.m, 2)))    #e_cache第一列给出e_cache是否有效的标志位；第二列给出实际E值
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:,i] = kernel_trans(self.data_mat, self.data_mat[i,:], k_tup)

def load_data_set(file_name):
    with open(file_name) as fd:
        fd_data = fd.readlines()

    data_set = []; label_set = []
    for data in fd_data:
        data = data.strip().split('\t')
        data = [float(value) for value in data]
        data_set.append(data[: -1])
        label_set.append(data[-1])
    return data_set, label_set

def select_j_rand(i, m):    #选择一个不等于i的j
    j = i
    while j == i:
        j = int(random.uniform(0, m))
    return j

def clip_alpha(aj, H, L):    #数值调整
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

"""计算E值并返回"""
def calc_Ek(o_S, k):    #传入参数：数据结构、索引
    """修改1"""
    f_Xk = float(multiply(o_S.alphas, o_S.label_mat).T * o_S.K[:,k] + o_S.b)
    #f_Xk = float(multiply(o_S.alphas, o_S.label_mat).T * (o_S.data_mat * o_S.data_mat[k,:].T)) + o_S.b
    Ek = f_Xk - float(o_S.label_mat[k])
    return Ek

"""选择第二个alpha值"""
def select_j(i, o_S, Ei):    #传入参数：索引i、数据结构、误差Ei
    max_K = -1; max_delta_E = 0; Ej = 0
    o_S.e_cache[i] = [1, Ei]    #将输入值Ei设置为“有效的”；有效的是指已经计算完毕
    valid_e_cache_list = nonzero(o_S.e_cache[:,0].A)[0]    #nonzero()返回e_cache[:,0].A中非零值的索引；.A表示将矩阵转化为数组
    if len(valid_e_cache_list) > 1:
        for k in valid_e_cache_list:
            if k == i:
                continue
            Ek = calc_Ek(o_S, k)
            delta_E = abs(Ei - Ek)
            if (delta_E > max_delta_E):
                max_K = k; max_delta_E = delta_E; Ej = Ek
        return max_K, Ej
    else:
        j = select_j_rand(i, o_S.m)
        Ej = calc_Ek(o_S, j)
    return j, Ej

"""更新误差值"""
def update_Ek(o_S, k):    #传入参数：数据结构、索引k
    Ek = calc_Ek(o_S, k)
    o_S.e_cache[k] = [1, Ek]

"""与简化版SMO几乎一样，除了使用数据结构作为输入、使用select_j更新j，且在alpha改变时更新e_cache"""
def inner_L(i, o_S):    #输入参数：索引、数据结构--与简化版SMO不同之处1
    Ei = calc_Ek(o_S, i)
    if ((o_S.label_mat[i] * Ei < -o_S.toler) and (o_S.alphas[i] < o_S.C)) or ((o_S.label_mat[i] * Ei > o_S.toler) and (o_S.alphas[i] > 0)):
        j, Ej = select_j(i, o_S, Ei)
        alpha_iold = o_S.alphas[i].copy();
        alpha_jold = o_S.alphas[j].copy();
        if (o_S.label_mat[i] != o_S.label_mat[j]):
            L = max(0, o_S.alphas[j] - o_S.alphas[i])
            H = min(o_S.C, o_S.C + o_S.alphas[j] - o_S.alphas[i])
        else:
            L = max(0, o_S.alphas[j] + o_S.alphas[i] - o_S.C)
            H = min(o_S.C, o_S.alphas[j] + o_S.alphas[i])
        if L == H:
            print("L==H")
            return 0
        """修改1"""
        eta = 2.0 * o_S.K[i, j] - o_S.K[i, i] - o_S.K[j, j]
        #eta = 2.0 * o_S.data_mat[i, :] * o_S.data_mat[j, :].T - o_S.data_mat[i, :] * o_S.data_mat[i, :].T - o_S.data_mat[j, :] * o_S.data_mat[j, :].T
        if (eta >= 0):
            print("eta >= 0")
            return 0
        o_S.alphas[j] -= o_S.label_mat[j] * (Ei - Ej) / eta
        o_S.alphas[j] = clip_alpha(o_S.alphas[j], H, L)
        update_Ek(o_S, j)    #与简化版SMO不同之处3
        if (abs(o_S.alphas[j] - alpha_jold) < 0.00001):
            print("j not moving enough")
            return 0
        o_S.alphas[i] += o_S.label_mat[j] * o_S.label_mat[i] * (alpha_jold - o_S.alphas[j])
        update_Ek(o_S, i)
        """修改2"""
        b1 = o_S.b - Ei - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.K[i, i] - \
            o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.K[i, j]
        b2 = o_S.b - Ej - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.K[i, i] - \
        o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.K[j, j]
        # b1 = o_S.b - Ei - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.data_mat[i, :] * o_S.data_mat[i, :].T -  \
        #      o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.data_mat[i, :] * o_S.data_mat[j, :].T
        # b2 = o_S.b - Ej - o_S.label_mat[i] * (o_S.alphas[i] - alpha_iold) * o_S.data_mat[i, :] * o_S.data_mat[j, :].T - \
        #      o_S.label_mat[j] * (o_S.alphas[j] - alpha_jold) * o_S.data_mat[j, :] * o_S.data_mat[j, :].T
        if (0 < o_S.alphas[i]) and (o_S.C > o_S.alphas[i]):
            o_S.b = b1
        elif (0 < o_S.alphas[j]) and (o_S.C > o_S.alphas[j]):
            o_S.b = b2
        else:
            o_S.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0

"""完整版SMO算法"""
def smop(data_set, label_set, C, toler, max_iter, k_tup):    #输入参数：数据集、标签集、常数C、容错率、最大循环次数
    o_S = optStruct(data_set, label_set, C, toler, k_tup)    #初始化数据结构
    iter = 0
    entire_set = True    #是否遍历完整个数据集
    alpha_pairs_changed = 0
    while (iter < max_iter) and ((alpha_pairs_changed > 0) or (entire_set)):
        alpha_pairs_changed = 0
        if entire_set:
            for i in range(o_S.m):
                alpha_pairs_changed += inner_L(i, o_S)
                print("full set, iter: %d i: %d, pairs changed %d" % (iter, i, alpha_pairs_changed))
        else:
            non_bound_is = nonzero((o_S.alphas.A > 0) * (o_S.alphas.A < C))[0]
            for i in non_bound_is:
                alpha_pairs_changed += inner_L(i, o_S)
                print("non-bound, iter: %d i: %d, pairs changed %d" % (iter, i, alpha_pairs_changed))
        iter += 1    #此处的迭代定义为一次循环过程
        if entire_set:
            entire_set = False
        elif (alpha_pairs_changed == 0):
            entire_set = True
        print("iteration number: %d" % iter)
    return o_S.b, o_S.alphas

"""计算W"""
def calc_ws(alphas, data_set, label_set):
    data_mat = mat(data_set)
    label_mat = mat(label_set).T
    m , n =shape(data_mat)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * label_mat[i], data_mat[i,:].T)
    return w

def kernel_trans(data_set, data_set_line, k_tup):
    data_mat = mat(data_set)
    m, n = shape(data_mat)
    K = mat(zeros((m, 1)))
    """此处只给出两种核函数的类型，可以继续添加"""
    if k_tup[0] == 'lin':    #线性核函数：内积计算在“所有数据集”和“数据集中的一行中展开”
        K = data_mat * data_set_line.T
    elif k_tup[0] == 'rbf':    #径向基核函数
        for j in range(m):
            delta_row = data_mat[j,:] - data_set_line
            K[j] = delta_row * delta_row.T
        K = exp(K / (-1 * k_tup[1]**2))    #numpy中/代表对元素展开即元素间的除法，而不像在MATLAB中一样计算矩阵的逆
    else:    #遇到无法识别的元组则退出
        raise NameError('Houston We Have a Problem -- That Kernel is not recogized')
    return K

def test_rbf(k1=1.3):
    data_set, label_set = load_data_set('C:/Users\Administrator\Python\pythonStudy\data/trainingDataSVM.txt')
    b, alphas = smop(data_set, label_set, 3, 0.0001, 10000, ('rbf', k1))
    data_mat = mat(data_set); label_mat = mat(label_set).T
    sv_ind = nonzero(alphas.A > 0)[0]
    data_sv = data_mat[sv_ind]; label_sv = label_mat[sv_ind]
    print("There are %d Support Vectors" % shape(data_sv)[0])
    m , n = shape(data_mat)
    error_count = 0
    for i in range(m):
        kernel_eval = kernel_trans(data_sv, data_mat[i,:], ('rbf', k1))
        predict = kernel_eval.T * multiply(label_sv, alphas[sv_ind]) + b
        if sign(predict) != sign(label_set[i]):
            error_count += 1
    print("The training error rate is: %f" % (float(error_count) / m))
    error_count = 0
    data_set, label_set = load_data_set('C:/Users\Administrator\Python\pythonStudy\data/testDataSVM.txt')
    data_mat = mat(data_set)
    m, n = shape(data_mat)
    for i in range(m):
        kernel_eval = kernel_trans(data_sv, data_mat[i,:], ('rbf', k1))
        predict = kernel_eval.T * multiply(label_sv, alphas[sv_ind]) + b
        if sign(predict) != sign(label_set[i]):
            error_count += 1
    print("The test error rate is: %f" % (float(error_count) / m))

if __name__ == '__main__':
    test_rbf(0.5)

---

1. $label\ast (W^{T}X+b)$称为点到分隔面的函数间隔，$label\ast (W^{T}X+b)\cdot \frac{1}{||W||}$称为点到分隔面的几何间隔。 ↩︎

2. John C.Platt, “Using Analytic QP and Sparseness to Speed Training of Support Vector Machines” in Advances in Neural Information Processing Systems 11, M.S.Kearns, S.A.Solla, D.A.Cohn, eds(MIT Press, 1999), 557-63 ↩︎