机器学习算法实战项目—支持向量机（1）—应用简化版SMO算法处理小规模数据集

叁歲伎倆 2023-10-02 20:24 10阅读 0赞

提供数据集如下:  
链接：https://pan.baidu.com/s/1\_k3vn7Rwt6DnKI\_VcO8MPw  
提取码：6666  
此外，完整版的学习笔记已经上传，有需要的可以去我的主页找，Latex排版后生成的文档。

## SVM的代码实现 ##

### 引言:Platt的SMO算法 ###

1996年，John Platt发布了一个称为SMO的强大算法，用于训练SVM。

SMO表示序列最小优化。Platt的SMO算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。在结果完全相同的同时，SMO算法的求解时间短很多。

SMO算法的目标是求出一系列 α \\alpha α和 b b b，一旦求出了这些 α \\alpha α，就很容易计算出权重向量 ω \\omega ω并得到分隔超平面。

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

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

首先在数据集上遍历每一个 α \\alpha α，然后在剩下的 α \\alpha α集合中随机选择另一个 α \\alpha α，从而构建 α \\alpha α对。这里有一点相当重要，就是我们要同时改变两个 α \\alpha α。之所以这样做是因为我们有一个约束条件：  
 ∑ i = 1 N α i y i = 0 \\sum\_\{i=1\}^\{N\}\\alpha\_iy\_i=0 i=1∑Nαiyi=0

此时我们使用的是线性支持向量机。

为此，我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们也需要另一个辅助函数，用于在数值太大时对其进行调整。下面给出了这两个函数的实现:

"""
    函数说明：读取数据,处理文本数据
    Parameters:
        file_name - 文件名
    Returns:
        data_mat - 数据矩阵
        label_mat - 数据标签
    """
    def load_dataset(file_name):
        # 数据矩阵
        data_mat = []
        # 标签向量
        label_mat = []
        # 打开文件
        fr = open(file_name)
        # 逐行读取
        for line in fr.readlines():
            # 去掉每一行首尾的空白符，例如'\n','\r','\t',' '
            # 将每一行内容根据'\t'符进行切片
            line_array = line.strip().split('\t')
            # 添加数据(100个元素排成一行)
            data_mat.append([float(line_array[0]), float(line_array[1])])
            # 添加标签(100个元素排成一行)
            label_mat.append(int(line_array[2]))
        return data_mat, label_mat

"""
    函数说明：随机选择alpha_j。alpha的选取，随机选择一个不等于i值的j
    Parameters:
        i - 第一个alpha的下标
        m - alpha参数个数 
    Returns:
        j - 返回选定的数字
    """
    def select_J_rand(i, m):
        """
        :param i: 表示alpha_i
        :param m: 表示样本数
        :return:
        """
        j = i
        while(j == i):
            # 如果i和j相等，那么就从样本中随机选取一个，可以认为j就是选择的alpha_2
            # uniform()方法将随机生成一个实数，它在[x, y)范围内
            j = int(random.uniform(0, m))
        return j

"""
    函数说明：修剪alpha
    Parameters:
        a_j - alpha值
        H - alpha上限
        L - alpha下限
    Returns:
        a_j - alpha值
    用于调整大于H或小于L的alpha值。
    """
    def clip_alpha(a_j, H, L):
        if a_j > H:
            a_j = H
        if L > a_j:
            a_j = L
        return a_j

第一个函数就是我们所熟知的load-dataset()函数，该函数打开文件并对其进行逐行解析，从而得到每行的类标签和整个数据矩阵。

下一个函数select-J-rand()有两个参数值，其中 i i i是第一个 α \\alpha α的下标， m m m是所有 α \\alpha α的数目。只要函数值不等于输入值 i i i，函数就会进行随机选择。

最后一个辅助函数就是clip-alpha()，它是用于调整大于 H H H或小于 L L L的 α \\alpha α值。

数据标签为:[-1, -1, 1, -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, 1,
    1......1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, -1, -1,
    -1, -1, -1]

可以看得出来，这里采用的类别标签是-1和1，而不是0和1。

我们可视化一下自己的数据集:

fig = plt.figure ()
    ax = fig.add_subplot(111)
    cm = mpl.colors.ListedColormap(['g', 'r'])
    ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm , s=20)

![51128e3026ca4e3c9ab6c8637146b5e1.png][]

我们的伪代码大致如下:

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

下面的代码是SMO算法的一个简单版本:

"""
    函数说明：简化版SMO算法
    Parameters:
        data_mat_in - 数据矩阵
        class_labels - 数据标签
        C - 松弛变量
        toler - 容错率
        max_iter - 最大迭代次数 
    Returns:
        None
    """
    def smo_simple(data_mat_in, class_labels, C, toler, max_iter):
        """
        :param data_mat_in: 相当于我们的x
        :param class_labels: 相当于我们的y
        :param C: 惩罚因子
        :param toler:
        :param max_iter:
        :return:
        """
        # 转换为numpy的mat矩阵存储(100,2)
        data_matrix = np.mat(data_mat_in) #相当于x
        # 转换为numpy的mat矩阵存储并转置(100,1)
        label_mat = np.mat(class_labels).transpose() #相当于y
        # 初始化b参数，统计data_matrix的维度,m:行；n:列
        b = 0
        # 统计dataMatrix的维度,m:100行；n:2列
        m, n = np.shape(data_matrix)
        # 初始化alpha参数，设为0
        alphas = np.mat(np.zeros((m, 1)))
        # 初始化迭代次数
        iter_num = 0
        # 最多迭代maxIter次
        while(iter_num < max_iter):
            alpha_pairs_changed = 0
            for i in range(m):
                # 步骤1：计算误差Ei
                # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
                # .T为转置，转置的目的是因为后面的核函数是一个向量。核函数我们最原始的方式，可以直接使用点乘就可以，即x_i.x
                fxi = float(np.multiply(alphas, label_mat).T * (data_matrix * data_matrix[i, :].T)) + b
                # 误差项计算公式
                Ei = fxi - float(label_mat[i])
                # 优化alpha，设定一定的容错率
                if((label_mat[i] * Ei < -toler) and (alphas[i] < C)) or ((label_mat[i] * Ei > toler) and (alphas[i] > 0)):
                    # 随机选择另一个alpha_i成对比优化的alpha_j
                    j = select_J_rand(i, m)
                    # 步骤1，计算误差Ej
                    fxj = float(np.multiply(alphas, label_mat).T * (data_matrix * data_matrix[j, :].T)) + b
                    # 误差项计算公式
                    Ej = fxj - float(label_mat[j])
                    # 保存更新前的alpha值，使用深拷贝(完全拷贝)A深层拷贝为B，A和B是两个独立的个体
                    alpha_i_old = alphas[i].copy()
                    alpha_j_old = alphas[j].copy()
                    # 步骤2：计算上下界H和L
                    if(label_mat[i] != label_mat[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
                    eta = 2.0 * data_matrix[i, :] * data_matrix[j, :].T - data_matrix[i, :] * data_matrix[i, :].T - data_matrix[j, :] * data_matrix[j, :].T
                    if eta >= 0:
                        print("eta>=0")
                        continue
                    # 步骤4：更新alpha_j
                    alphas[j] -= label_mat[j] * (Ei - Ej) / eta
                    # 步骤5：修剪alpha_j
                    alphas[j] = clip_alpha(alphas[j], H, L)
                    if(abs(alphas[j] - alpha_j_old) < 0.00001):
                        print("alpha_j变化太小")
                        continue
                    # 步骤6：更新alpha_i
                    alphas[i] += label_mat[j] * label_mat[i] * (alpha_j_old - alphas[j])
                    # 步骤7：更新b_1和b_2
                    b1 = b - Ei - label_mat[i] * (alphas[i] - alpha_i_old) * data_matrix[i, :] * data_matrix[i, :].T - label_mat[j] * (alphas[j] - alpha_j_old) * data_matrix[j, :] * data_matrix[i, :].T
                    b2 = b - Ej - label_mat[i] * (alphas[i] - alpha_i_old) * data_matrix[i, :] * data_matrix[j, :].T - label_mat[j] * (alphas[j] - alpha_j_old) * data_matrix[j, :] * data_matrix[j, :].T
                    # 步骤8：根据b_1和b_2更新b
                    if(0 < alphas[i] < C):
                        b = b1
                    elif(0 < alphas[j] < C):
                        b = b2
                    else:
                        b = (b1 + b2) / 2.0
                    # 统计优化次数
                    alpha_pairs_changed += 1
                    # 打印统计信息
                    print("第%d次迭代 样本：%d， alpha优化次数：%d" % (iter_num, i, alpha_pairs_changed))
            # 更新迭代次数
            if(alpha_pairs_changed == 0):
                iter_num += 1
            else:
                iter_num = 0
            print("迭代次数：%d" % iter_num)
        return b, alphas

该函数有5个输入参数，分别是：数据集、类别标签、常数C、容错率和退出前最大的循环次数。

上述函数将多个列表和输入参数转换成NumPy矩阵，这样就可以简化很多数学处理操作。由于转置了类别标签，因此我们得到的就是一个列向量而不是列表。于是类别标签向量的每行元素都和数据矩阵中的行一一对应。

我们也可以通过矩阵data-mat-in的shape属性得到常数m和n。最后，我们就可以构建一个alpha列矩阵，矩阵中元素都初始化为0，并建立一个iter变量。该变量存储的则是在没有任何alpha改变  
的情况下遍历数据集的次数。当该变量达到输入值max-iter时，函数结束运行并退出。

每次循环当中，将alpha-pairs-changed先设为0，然后再对整个集合顺序遍历。变量alpha-pairs-changed用于记录alpha是否已经进行优化。当然，在循环结束时就会得知这一点。

首先，fXi能够计算出来，这就是我们预测的类别。然后，基于这个实例的预测结果和真实结果的比对，就可以计算误差Ei。如果误差很大，那么可以对该数据实例所对应的alpha值进行优化。

在if语句中，不管是正间隔还是负间隔都会被测试。并且在该if语句中，也要同时检查alpha值，以保证其不能等于0或C。由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值的话，那么它们就已经在“边界”上了，因而不再能够减小或增大，因此也就不值得再对它们进行优化了。

接下来，可以利用程序中的辅助函数来随机选择第二个alpha值，即alpha\[j\] 。

同样，可以采用第一个alpha值（alpha\[i\]）的误差计算方法，来计算这个alpha值的误差。

这个过程可以通过copy()的方法来实现，因此稍后可以将新的alpha值与老的alpha值进行比较。Python则会通过引用的方式传递所有列表，所以必须明确地告知Python要为alpha-i-old和alpha-j-old分配新的内存；否则的话，在对新值和旧值进行比较时，我们就看不到新旧值的变化。

之后我们开始计算L和H，它们用于将alpha\[j\]调整到0到C之间。如果L和H相等，就不做任何改变，直接执行continue语句。这在Python中，则意味着本次循环结束直接运行下一次for的循环。

Eta是alpha\[j\]的最优修改量，在那个很长的计算代码行中得到。如果eta为0，那就是说需要退出for循环的当前迭代过程。该过程对真实SMO算法进行了简化处理。如果eta为0，那么计算新的alpha\[j\]就比较麻烦了，这里我们就不对此进行详细的介绍了。有需要的可以阅读Platt的原文来了解更多的细节。现实中，这种情况并不常发生，因此忽略这一部分通常也无伤大雅。

于是，可以计算出一个新的alpha\[j\]，然后利用程序中的辅助函数以及L与H值对其进行调整。

然后，就是需要检查alpha\[j\]是否有轻微改变。如果是的话，就退出for循环。然后，alpha\[i\]和alpha\[j\]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反（即如果一个增加，那么另外一个减少）。在对alpha\[i\]和alpha\[j\]进行优化之后，给这两个alpha值设置一个常数项b  
。

最后，在优化过程结束的同时，必须确保在合适的时机结束循环。如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha，同时可以增加alpha-pairs-changed的值。在for循环之外，需要检查alpha值是否做了更新，如果有更新则将iter设为0后继续运行程序。只有在所有数据集上遍历max-iter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环。

我们得到的计算 ω \\omega ω的程序如下:

"""
    函数说明：计算w
    Returns:
        data_mat - 数据矩阵
        label_mat - 数据标签
        alphas - alphas值
    Returns:
        w - 直线法向量
    """
    def get_w(data_mat_in, class_labels, alphas):
        """
        :param data_mat: x
        :param label_mat: y
        :param alphas:
        :return:
        """
        data_mat=np.mat(data_mat_in)
        label_mat=np.mat(class_labels).transpose()
        m,n=np.shape(data_mat)
        # 初始化w都为1
        w=np.zeros((n,1))
        # dot()函数是矩阵乘，而*则表示逐个元素相乘
        # w = sum(alpha_i * yi * xi)
        #循环计算
        for i in range(m):
            w+=np.multiply(alphas[i]*label_mat[i],data_mat[i,:].T)
        return w

我们得到参数的值分别为:

b= [[-3.83678533]]
    w= [[ 0.81412243]
     [-0.27106553]]
    alphas= [[0.12694755 0.24149422 0.36844176]]

我们的迭代过程如下所示:

.......此处省略部分结果
    迭代次数：5
    alpha_j变化太小
    迭代次数：6
    alpha_j变化太小
    迭代次数：7
    alpha_j变化太小
    迭代次数：8
    alpha_j变化太小
    迭代次数：9
    alpha_j变化太小
    迭代次数：10
    alpha_j变化太小
    迭代次数：11
    alpha_j变化太小
    迭代次数：12
    alpha_j变化太小
    迭代次数：13
    alpha_j变化太小
    迭代次数：14
    alpha_j变化太小
    迭代次数：15
    alpha_j变化太小
    迭代次数：16
    alpha_j变化太小
    迭代次数：17
    alpha_j变化太小
    迭代次数：18
    alpha_j变化太小
    迭代次数：19
    alpha_j变化太小
    迭代次数：20
    alpha_j变化太小
    迭代次数：21
    alpha_j变化太小
    迭代次数：22
    alpha_j变化太小
    迭代次数：23
    alpha_j变化太小
    迭代次数：24
    alpha_j变化太小
    迭代次数：25
    alpha_j变化太小
    迭代次数：26
    alpha_j变化太小
    迭代次数：27
    alpha_j变化太小
    迭代次数：28
    alpha_j变化太小
    迭代次数：29
    alpha_j变化太小
    迭代次数：30
    alpha_j变化太小
    迭代次数：31
    alpha_j变化太小
    迭代次数：32
    alpha_j变化太小
    迭代次数：33
    alpha_j变化太小
    迭代次数：34
    alpha_j变化太小
    迭代次数：35
    alpha_j变化太小
    迭代次数：36
    alpha_j变化太小
    迭代次数：37
    alpha_j变化太小
    迭代次数：38
    alpha_j变化太小
    迭代次数：39
    alpha_j变化太小
    迭代次数：40

为了解哪些数据点是支持向量，输入：

for i in range(100):
         if alphas[i]>0.0:
            print(data_mat[i],label_mat[i])

得到的结果如下：

[4.658191, 3.507396] -1
    [3.457096, -0.082216] -1
    [6.080573, 0.418886] 1

因为是线性可分，我们画出决策边界:

x = np.arange(-2.0, 12, 0.1)
    # w0x0+w1x1=0,x1=-w0x0/w1
    y = (-w[0] * x - b) / w[1]
    ax.plot(x,y.reshape(-1,1))
    ax.axis([-2,12,-8.6,7])
    plt.show()

![575acd56f7094d29afd79032f61ad5bd.png][]

我们希望同时找出支持向量，下面是书中记录的案例:

"""
    函数说明：分类结果可视化
    Returns:
        dataMat - 数据矩阵
        w - 直线法向量
        b - 直线截距
    Returns:
        None
    """
    def show_classifer(data_mat, w, b):
        data_mat, label_mat = load_dataset('testSet.txt')
        fig = plt.figure()
        ax = fig.add_subplot(111)
        cm = mpl.colors.ListedColormap(['g', 'r'])
        ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm, s=20)
        # 绘制直线
        x1 = max(data_mat)[0]
        x2 = min(data_mat)[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])
        # 找出支持向量点
        # enumerate在字典上是枚举、列举的意思
        for i, alpha in enumerate(alphas):
            # 支持向量机的点
            if(abs(alpha) > 0):
                x, y = data_mat[i]
                plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolors='red')
        plt.show()

![e19bb65bf2cc433c8001b9996f12a78e.png][]

在更大的数据集上，收敛时间会变得更长。在下一节中，我们优化算法，以提高算法的运行效率。

### 完整版代码如下： ###

# -*- coding: utf-8 -*-
    """
    Created on Sun Jul 18 2022
    @author: wzk
    """
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import numpy as np
    import random
    
    """
    函数说明：读取数据,处理文本数据
    Parameters:
        file_name - 文件名
    Returns:
        data_mat - 数据矩阵
        label_mat - 数据标签
    """
    def load_dataset(file_name):
        # 数据矩阵
        data_mat = []
        # 标签向量
        label_mat = []
        # 打开文件
        fr = open(file_name)
        # 逐行读取
        for line in fr.readlines():
            # 去掉每一行首尾的空白符，例如'\n','\r','\t',' '
            # 将每一行内容根据'\t'符进行切片
            line_array = line.strip().split('\t')
            # 添加数据(100个元素排成一行)
            data_mat.append([float(line_array[0]), float(line_array[1])])
            # 添加标签(100个元素排成一行)
            label_mat.append(int(line_array[2]))
        return data_mat, label_mat
    
    
    """
    函数说明：随机选择alpha_j。alpha的选取，随机选择一个不等于i值的j
    Parameters:
        i - 第一个alpha的下标
        m - alpha参数个数 
    Returns:
        j - 返回选定的数字
    """
    def select_J_rand(i, m):
        """
        :param i: 表示alpha_i
        :param m: 表示样本数
        :return:
        """
        j = i
        while(j == i):
            # 如果i和j相等，那么就从样本中随机选取一个，可以认为j就是选择的alpha_2
            # uniform()方法将随机生成一个实数，它在[x, y)范围内
            j = int(random.uniform(0, m))
        return j
    
    """
    函数说明：修剪alpha
    Parameters:
        a_j - alpha值
        H - alpha上限
        L - alpha下限
    Returns:
        a_j - alpha值
    用于调整大于H或小于L的alpha值。
    """
    def clip_alpha(a_j, H, L):
        if a_j > H:
            a_j = H
        if L > a_j:
            a_j = L
        return a_j
    
    """
    函数说明：简化版SMO算法
    Parameters:
        data_mat_in - 数据矩阵
        class_labels - 数据标签
        C - 松弛变量
        toler - 容错率
        max_iter - 最大迭代次数 
    Returns:
        None
    """
    def smo_simple(data_mat_in, class_labels, C, toler, max_iter):
        """
        :param data_mat_in: 相当于我们的x
        :param class_labels: 相当于我们的y
        :param C: 惩罚因子
        :param toler:
        :param max_iter:
        :return:
        """
        # 转换为numpy的mat矩阵存储(100,2)
        data_matrix = np.mat(data_mat_in) #相当于x
        # 转换为numpy的mat矩阵存储并转置(100,1)
        label_mat = np.mat(class_labels).transpose() #相当于y
        # 初始化b参数，统计data_matrix的维度,m:行；n:列
        b = 0
        # 统计dataMatrix的维度,m:100行；n:2列
        m, n = np.shape(data_matrix)
        # 初始化alpha参数，设为0
        alphas = np.mat(np.zeros((m, 1)))
        # 初始化迭代次数
        iter_num = 0
        # 最多迭代maxIter次
        while(iter_num < max_iter):
            alpha_pairs_changed = 0
            for i in range(m):
                # 步骤1：计算误差Ei
                # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
                # .T为转置，转置的目的是因为后面的核函数是一个向量。核函数我们最原始的方式，可以直接使用点乘就可以，即x_i.x
                fxi = float(np.multiply(alphas, label_mat).T * (data_matrix * data_matrix[i, :].T)) + b
                # 误差项计算公式
                Ei = fxi - float(label_mat[i])
                # 优化alpha，设定一定的容错率
                if((label_mat[i] * Ei < -toler) and (alphas[i] < C)) or ((label_mat[i] * Ei > toler) and (alphas[i] > 0)):
                    # 随机选择另一个alpha_i成对比优化的alpha_j
                    j = select_J_rand(i, m)
                    # 步骤1，计算误差Ej
                    fxj = float(np.multiply(alphas, label_mat).T * (data_matrix * data_matrix[j, :].T)) + b
                    # 误差项计算公式
                    Ej = fxj - float(label_mat[j])
                    # 保存更新前的alpha值，使用深拷贝(完全拷贝)A深层拷贝为B，A和B是两个独立的个体
                    alpha_i_old = alphas[i].copy()
                    alpha_j_old = alphas[j].copy()
                    # 步骤2：计算上下界H和L
                    if(label_mat[i] != label_mat[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
                    eta = 2.0 * data_matrix[i, :] * data_matrix[j, :].T - data_matrix[i, :] * data_matrix[i, :].T - data_matrix[j, :] * data_matrix[j, :].T
                    if eta >= 0:
                        print("eta>=0")
                        continue
                    # 步骤4：更新alpha_j
                    alphas[j] -= label_mat[j] * (Ei - Ej) / eta
                    # 步骤5：修剪alpha_j
                    alphas[j] = clip_alpha(alphas[j], H, L)
                    if(abs(alphas[j] - alpha_j_old) < 0.00001):
                        print("alpha_j变化太小")
                        continue
                    # 步骤6：更新alpha_i
                    alphas[i] += label_mat[j] * label_mat[i] * (alpha_j_old - alphas[j])
                    # 步骤7：更新b_1和b_2
                    b1 = b - Ei - label_mat[i] * (alphas[i] - alpha_i_old) * data_matrix[i, :] * data_matrix[i, :].T - label_mat[j] * (alphas[j] - alpha_j_old) * data_matrix[j, :] * data_matrix[i, :].T
                    b2 = b - Ej - label_mat[i] * (alphas[i] - alpha_i_old) * data_matrix[i, :] * data_matrix[j, :].T - label_mat[j] * (alphas[j] - alpha_j_old) * data_matrix[j, :] * data_matrix[j, :].T
                    # 步骤8：根据b_1和b_2更新b
                    if(0 < alphas[i] < C):
                        b = b1
                    elif(0 < alphas[j] < C):
                        b = b2
                    else:
                        b = (b1 + b2) / 2.0
                    # 统计优化次数
                    alpha_pairs_changed += 1
                    # 打印统计信息
                    print("第%d次迭代 样本：%d， alpha优化次数：%d" % (iter_num, i, alpha_pairs_changed))
            # 更新迭代次数
            if(alpha_pairs_changed == 0):
                iter_num += 1
            else:
                iter_num = 0
            print("迭代次数：%d" % iter_num)
        return b, alphas
    
    """
    函数说明：计算w
    Returns:
        data_mat - 数据矩阵
        label_mat - 数据标签
        alphas - alphas值
    Returns:
        w - 直线法向量
    """
    def get_w(data_mat_in, class_labels, alphas):
        """
        :param data_mat: x
        :param label_mat: y
        :param alphas:
        :return:
        """
        data_mat=np.mat(data_mat_in)
        label_mat=np.mat(class_labels).transpose()
        m,n=np.shape(data_mat)
        # 初始化w都为1
        w=np.zeros((n,1))
        # dot()函数是矩阵乘，而*则表示逐个元素相乘
        # w = sum(alpha_i * yi * xi)
        #循环计算
        for i in range(m):
            w+=np.multiply(alphas[i]*label_mat[i],data_mat[i,:].T)
        return w
    
    """
    函数说明：分类结果可视化
    Returns:
        dataMat - 数据矩阵
        w - 直线法向量
        b - 直线截距
    Returns:
        None
    """
    def show_classifer(data_mat, w, b):
        data_mat, label_mat = load_dataset('testSet.txt')
        fig = plt.figure()
        ax = fig.add_subplot(111)
        cm = mpl.colors.ListedColormap(['g', 'r'])
        ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm, s=20)
        # 绘制直线
        x1 = max(data_mat)[0]
        x2 = min(data_mat)[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])
        # 找出支持向量点
        # enumerate在字典上是枚举、列举的意思
        for i, alpha in enumerate(alphas):
            # 支持向量机的点
            if(abs(alpha) > 0):
                x, y = data_mat[i]
                plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolors='red')
        plt.show()
    
    if __name__ == '__main__':
        data_mat, label_mat = load_dataset('testSet.txt')
        print('数据标签为:\n',label_mat)
        fig = plt.figure()
        ax = fig.add_subplot(111)
        cm = mpl.colors.ListedColormap(['g', 'r'])
        ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm, s=20)
    
        b, alphas = smo_simple(data_mat, label_mat, 0.6, 0.001, 40)
        w = get_w(data_mat, label_mat, alphas)
        print('b=',b)
        print('w=',w)
        print('alphas=',alphas[alphas>0])
        for i in range(100):
            if alphas[i]>0.0:
                print(data_mat[i],label_mat[i])
        x = np.arange(-2.0, 12, 0.1)
        # w0x0+w1x1=0,x1=-w0x0/w1
        y = (-w[0] * x - b) / w[1]
        ax.plot(x,y.reshape(-1,1))
        ax.axis([-2,12,-8.6,7])
        plt.show()
        show_classifer(data_mat,w,b)

[51128e3026ca4e3c9ab6c8637146b5e1.png]: https://img-blog.csdnimg.cn/51128e3026ca4e3c9ab6c8637146b5e1.png
[575acd56f7094d29afd79032f61ad5bd.png]: https://img-blog.csdnimg.cn/575acd56f7094d29afd79032f61ad5bd.png
[e19bb65bf2cc433c8001b9996f12a78e.png]: https://img-blog.csdnimg.cn/e19bb65bf2cc433c8001b9996f12a78e.png