机器学习算法实战项目—支持向量机（2）—完整版的SMO算法

小灰灰 2023-09-28 16:19 4阅读 0赞

### 2.完整版的SMO算法 ###

在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有什么问题的，但是在更大的数据集上的运行速度就会变慢。

刚才已经讨论了简化版SMO算法，下面我们就讨论完整版的Platt SMO算法。

在这两个版本中，实现alpha的更改和代数运算的优化环节一模一样。在优化过程中，唯一的不同就是选择alpha的方式。

完整版的Platt SMO算法应用了一些能够提速的启发方法。上一节的例子在执行时存在一定的时间提升空间。

Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。

而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。

对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。

在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。

在简化版SMO算法中，我们会在选择j之后计算错误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。

辅助函数和上面的步骤一样，即下面的三个函数:

"""
    函数说明：读取数据
    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(float(line_array[2]))
        return data_mat, label_mat
    
    """
    函数说明：随机选择alpha_j
    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    def select_j_random(i, m):
        j = i
        while(j == i):
            # uniform()方法将随机生成一个实数，它在[x, y)范围内
            j = int(random.uniform(0, m))
        return j
    
    """
    函数说明：修剪alpha_j
    Parameters:
        aj - alpha_j值
        H - alpha上限
        L - alpha下限 
    Returns:
        aj - alpha_j值
    """
    def clip_alpha(aj, H, L):
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj

我们首先定义一个数据结构:

"""
    类说明：维护所有需要操作的值
    Parameters:
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率   
    Returns:
        None
    """
    """
    定义一个新的数据结构
    """
    class opt_struct:
        def __init__(self, data_mat_in, class_labels, C, toler):
            # 数据矩阵
            self.X = data_mat_in #传进来的数据
            # 数据标签
            self.label_mat = class_labels
            # 松弛变量
            self.C = C
            # 容错率
            self.tol = toler
            # 矩阵的行数
            self.m = np.shape(data_mat_in)[0]
            # 根据矩阵行数初始化alphas矩阵，一个m行1列的全零列向量
            self.alphas = np.mat(np.zeros((self.m, 1)))
            # 初始化b参数为0
            self.b = 0
            # 根据矩阵行数初始化误差缓存矩阵，第一列为是否有效标志位，其中0无效，1有效;第二列为实际的误差Ei的值
            """
            我们之前的定义为:Ei=gxi-yi
            yi是标签的实际值。
            gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
            因为误差值经常用到，所以希望每次计算后放到一个缓存当中，将ecache一分为二，第一列是标志位，取值为0或者1，为1时表示已经算出来
            """
            self.ecache = np.mat(np.zeros((self.m, 2)))

首要的事情就是建立一个数据结构来保存所有的重要值，而这个过程可以通过一个对象来完成。只是作为一个数据结构来使用对象。

在将值传给函数时，我们可以通过将所有数据移到一个结构中来实现，这样就可以省掉手工输入的麻烦了。而此时，数据就可以通过一个对象来进行传递。实际上，当完成其实现时，可以很容易通过Python的字典来完成。

该方法可以实现其成员变量的填充。除了增加了一个m×2的矩阵成  
员变量ecache之外，这些做法和简化版SMO一模一样。ecache的第一列给出的是ecache是否有效的标志位，而第二列给出的是实际的E值。

对于给定的alpha值，第一个辅助函数cal-Ek()能够计算E值并返回。以前，该过程是采用内嵌的方式来完成的，但是由于该过程在这个版本的SMO算法中出现频繁，这里必须要将其单独拎出来。

"""
    函数说明：计算误差
    Parameters:
        os - 数据结构
        k - 标号为k的数据
    Returns:
        Ek - 标号为k的数据误差
    """
    def cal_Ek(os, k):
        # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
        # .T为转置
        fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
        # 计算误差项
        Ek = fXk - float(os.label_mat[k])
        # 返回误差项
        return Ek

函数select-j()用于选择第二个alpha或者说内循环的alpha值。回想一下，这里的目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关。

首先将输入值Ei在缓存中设置成为有效的。这里的有效意味着它已经计算好了。在ecache中，代码nonzero(os.ecache\[:,0\].A)\[0\]构建出了一个非零表。

NumPy函数nonzero()返回了一个列表，而这个列表中包含以输入列表为目录的列表值，当然这里的值并非零。nonzero()语句返回的是非零E值所对应的alpha值，而不是E值本身。

程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果这是第一次循环的话，那么就随机选择一个alpha值。当然也存在有许多更复杂的方式来处理第一次循环的情况，而上述做法就能够满足我们的目的。

"""
    函数说明：内循环启发方式2
    选择第二个待优化的alpha_j，选择一个误差最大的alpha_j
    即，我们在选择alpha_2的时候做了改进，选择误差最大的
    Parameters:
        i - 标号为i的数据的索引值
        oS - 数据结构
        Ei - 标号为i的数据误差 
    Returns:
        j - 标号为j的数据的索引值
        maxK - 标号为maxK的数据的索引值
        Ej - 标号为j的数据误差
    """
    def select_j(i, os, Ei):
        # 初始化
        max_K = -1 #下标的索引值
        max_delta_E = 0
        Ej = 0
        # 根据Ei更新误差缓存，即先计算alpha_1以及E1值
        os.ecache[i] = [1, Ei] #放入缓存当中，设为有效
        # 对一个矩阵.A转换为Array类型
        # 返回误差不为0的数据的索引值
        valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
        # 有不为0的误差
        if(len(valid_ecache_list) > 1):
            # 遍历，找到最大的Ek
            for k in valid_ecache_list: #迭代所有有效的缓存，找到误差最大的E
                # 不计算k==i节省时间
                if k == i: #不选择和i相等的值
                    continue
                # 计算Ek
                Ek = cal_Ek(os, k)
                # 计算|Ei - Ek|
                delta_E = abs(Ei - Ek)
                # 找到maxDeltaE
                if(delta_E > max_delta_E):
                    max_K = k
                    max_delta_E = delta_E
                    Ej = Ek
            # 返回max_K，Ej
            return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
        # 没有不为0的误差
        else: #第一次循环时是没有有效的缓存值的，所以随机选一个（仅会执行一次）
            # 随机选择alpha_j的索引值
            j = select_j_random(i, os.m)
            # 计算Ej
            Ej = cal_Ek(os, j)
        # 返回j，Ej
        return j, Ej

最后一个辅助函数是update-Ek()，它会计算误差值并存入缓存当中。在对alpha值进行优化之后会用到这个值。

"""
    函数说明：计算Ek,并更新误差缓存
    Parameters:
        os - 数据结构
        k - 标号为k的数据的索引值
    Returns:
        None
    """
    def update_Ek(os, k):
        # 计算Ek
        Ek = cal_Ek(os, k)
        # 更新误差缓存
        os.ecache[k] = [1, Ek]

接下来将简单介绍一下用于寻找决策边界的优化例程:

"""
    函数说明：优化的SMO算法
    Parameters:
        i - 标号为i的数据的索引值
        os - 数据结构
    Returns:
        1 - 有任意一对alpha值发生变化
        0 - 没有任意一对alpha值发生变化或变化太小
    """
    def innerL(i, os):
        # 步骤1：计算误差Ei
        Ei = cal_Ek(os, i)
        # 优化alpha,设定一定的容错率
        if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
            # 使用内循环启发方式2选择alpha_j,并计算Ej
            j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
            # 保存更新前的alpha值，使用深层拷贝
            alpha_i_old = os.alphas[i].copy()
            alpha_j_old = os.alphas[j].copy()
            # 步骤2：计算上界H和下界L
            if(os.label_mat[i] != os.label_mat[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.label_mat[j] * (Ei - Ej) / eta
            # 步骤5：修剪alpha_j
            os.alphas[j] = clip_alpha(os.alphas[j], H, L)
            # 更新Ej至误差缓存
            update_Ek(os, j)
            if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
                print("alpha_j变化太小")
                return 0
            # 步骤6：更新alpha_i
            os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
            # 更新Ei至误差缓存
            update_Ek(os, i)
            # 步骤7：更新b_1和b_2:
            b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
            b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
            # 步骤8：根据b_1和b_2更新b
            if(0 < os.alphas[i] < os.C):
                os.b = b1
            elif(0 < os.alphas[j] < os.C):
                os.b = b2
            else:
                os.b = (b1 + b2) / 2.0
            return 1 #表示有更新
        else:
            return 0 #表示没有更新

代码几乎与上一节中给出的smo-simple函数一摸一样，但是这里的代码已经使用了自己的数据结构。该结构在参数os中传递。

完整版SMO算法的外循环代码:

"""
    函数说明：完整的线性SMO算法
    Parameters:
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
        maxIter - 最大迭代次数
    Returns:
        oS.b - SMO算法计算的b
        oS.alphas - SMO算法计算的alphas
    """
    def smo_p(data_mat_in, class_labels, C, toler, max_iter):
        # 初始化数据结构
        os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
        # 初始化当前迭代次数
        iter = 0
        entrie_set = True #是否在全部数据集上迭代
        alpha_pairs_changed = 0
        # 遍历整个数据集alpha都没有更新或者超过最大迭代次数，则退出循环
        while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
            alpha_pairs_changed = 0
            if entrie_set: # 遍历整个数据集
                for i in range(os.m):
                    # 使用优化的SMO算法
                    alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
                    print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
                iter += 1
            # 遍历非边界值
            else:
                # 遍历不在边界0和C的alpha
                non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
                for i in non_bound_i_s:
                    alpha_pairs_changed += innerL(i, os)
                    print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
                iter += 1
            # 遍历一次后改为非边界遍历
            if entrie_set:
                entrie_set = False #进行切换，遍历非边界数据集
            # 如果alpha没有更新，计算全样本遍历
            elif(alpha_pairs_changed == 0):
                entrie_set = True
            print("迭代次数:%d" % iter)
        # 返回SMO算法计算的b和alphas
        return os.b, os.alphas

其输入和函数smo-simple()完全一样。函数一开始构建一个数据结构来容纳所有的数据，然后需要对控制函数退出的一些变量进行初始化。整个代码的主体是while循环，这与smo-simple()有些类似，但是这里的循环退出条件更多一些。

当迭代次数超过指定的最大值，或者遍历整个集合都未对任意alpha对进行修改时，就退出循环。

这里的max-iter变量和函数smo-simple()中的作用有一点不同，后者当没有任何alpha发生改变时会将整个集合的一次遍历过程计成一次迭代，而这里的一次迭代定义为一次循环过程，而不管该循环具体做了什么事。此时，如果在优化过程中存在波动就会停止，因此这里的做法优于smo-simple()函数中的计数方法。

while循环的内部与smo-simple()中有所不同，一开始的for循环在数据集上遍历任意可能的alpha。

我们通过调用innerL()来选择第二个alpha，并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么会返回  
1。第二个for循环遍历所有的非边界alpha值，也就是不在边界0或C上的值。

接下来，我们对for循环在非边界循环和完整遍历之间进行切换，并打印出迭代次数。最后程序将会返回常数b和alpha值。

我们的部分结果如下:

全样本遍历:第0次迭代 样本:0, alpha优化次数:1
    全样本遍历:第0次迭代 样本:1, alpha优化次数:1
    全样本遍历:第0次迭代 样本:2, alpha优化次数:1
    全样本遍历:第0次迭代 样本:3, alpha优化次数:1
    全样本遍历:第0次迭代 样本:4, alpha优化次数:1
    全样本遍历:第0次迭代 样本:5, alpha优化次数:1
    全样本遍历:第0次迭代 样本:6, alpha优化次数:1
    全样本遍历:第0次迭代 样本:7, alpha优化次数:1
    全样本遍历:第0次迭代 样本:8, alpha优化次数:2
    全样本遍历:第0次迭代 样本:9, alpha优化次数:2
    ......
    非边界遍历:第1次迭代 样本:26, alpha优化次数:3
    alpha_j变化太小
    非边界遍历:第1次迭代 样本:29, alpha优化次数:3
    alpha_j变化太小
    非边界遍历:第1次迭代 样本:46, alpha优化次数:3
    alpha_j变化太小
    非边界遍历:第1次迭代 样本:54, alpha优化次数:3
    非边界遍历:第1次迭代 样本:55, alpha优化次数:3
    迭代次数:2
    alpha_j变化太小
    非边界遍历:第2次迭代 样本:8, alpha优化次数:0
    alpha_j变化太小
    非边界遍历:第2次迭代 样本:10, alpha优化次数:0
    非边界遍历:第2次迭代 样本:23, alpha优化次数:0
    alpha_j变化太小
    非边界遍历:第2次迭代 样本:29, alpha优化次数:0
    alpha_j变化太小
    非边界遍历:第2次迭代 样本:46, alpha优化次数:0
    alpha_j变化太小
    非边界遍历:第2次迭代 样本:54, alpha优化次数:0
    非边界遍历:第2次迭代 样本:55, alpha优化次数:0
    ......
    全样本遍历:第3次迭代 样本:97, alpha优化次数:0
    全样本遍历:第3次迭代 样本:98, alpha优化次数:0
    全样本遍历:第3次迭代 样本:99, alpha优化次数:0
    迭代次数:4

我们最终的决策边界和支持向量可视化为:  
![f03dd12888eb46f6bf1d4dd60dd0dfb5.png][]

b= [[-3.06127512]]
    w= [[ 0.69141217]
     [-0.3411598 ]]
    alphas= [[0.06158536 0.08744873 0.09509032 0.07285032 0.01603246 0.15524162]]
    [3.542485, 1.977398] -1.0
    [3.125951, 0.293251] -1.0
    [4.658191, 3.507396] -1.0
    [8.197181, 1.545132] 1.0
    [5.286862, -2.358286] 1.0
    [6.080573, 0.418886] 1.0

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

# -*- coding: utf-8 -*-
    import matplotlib.pyplot as plt
    import numpy as np
    import random
    import matplotlib as mpl
    """
    改进SMO算法以加快我们的SVM运行速度！
    """
    
    """
    函数说明：读取数据
    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(float(line_array[2]))
        return data_mat, label_mat
    
    """
    函数说明：随机选择alpha_j
    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    def select_j_random(i, m):
        j = i
        while(j == i):
            # uniform()方法将随机生成一个实数，它在[x, y)范围内
            j = int(random.uniform(0, m))
        return j
    
    """
    函数说明：修剪alpha_j
    Parameters:
        aj - alpha_j值
        H - alpha上限
        L - alpha下限 
    Returns:
        aj - alpha_j值
    """
    def clip_alpha(aj, H, L):
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj
    
    """
    类说明：维护所有需要操作的值
    Parameters:
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率   
    Returns:
        None
    """
    """
    定义一个新的数据结构
    """
    class opt_struct:
        def __init__(self, data_mat_in, class_labels, C, toler):
            # 数据矩阵
            self.X = data_mat_in #传进来的数据
            # 数据标签
            self.label_mat = class_labels
            # 松弛变量
            self.C = C
            # 容错率
            self.tol = toler
            # 矩阵的行数
            self.m = np.shape(data_mat_in)[0]
            # 根据矩阵行数初始化alphas矩阵，一个m行1列的全零列向量
            self.alphas = np.mat(np.zeros((self.m, 1)))
            # 初始化b参数为0
            self.b = 0
            # 根据矩阵行数初始化误差缓存矩阵，第一列为是否有效标志位，其中0无效，1有效;第二列为实际的误差Ei的值
            """
            我们之前的定义为:Ei=gxi-yi
            yi是标签的实际值。
            gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
            因为误差值经常用到，所以希望每次计算后放到一个缓存当中，将ecache一分为二，第一列是标志位，取值为0或者1，为1时表示已经算出来
            """
            self.ecache = np.mat(np.zeros((self.m, 2)))
    
    """
    函数说明：计算误差
    Parameters:
        os - 数据结构
        k - 标号为k的数据
    Returns:
        Ek - 标号为k的数据误差
    """
    def cal_Ek(os, k):
        # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
        # .T为转置
        fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
        # 计算误差项
        Ek = fXk - float(os.label_mat[k])
        # 返回误差项
        return Ek
    
    
    """
    函数说明：内循环启发方式2
    选择第二个待优化的alpha_j，选择一个误差最大的alpha_j
    即，我们在选择alpha_2的时候做了改进，选择误差最大的
    Parameters:
        i - 标号为i的数据的索引值
        oS - 数据结构
        Ei - 标号为i的数据误差 
    Returns:
        j - 标号为j的数据的索引值
        maxK - 标号为maxK的数据的索引值
        Ej - 标号为j的数据误差
    """
    def select_j(i, os, Ei):
        # 初始化
        max_K = -1 #下标的索引值
        max_delta_E = 0
        Ej = 0
        # 根据Ei更新误差缓存，即先计算alpha_1以及E1值
        os.ecache[i] = [1, Ei] #放入缓存当中，设为有效
        # 对一个矩阵.A转换为Array类型
        # 返回误差不为0的数据的索引值
        valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
        # 有不为0的误差
        if(len(valid_ecache_list) > 1):
            # 遍历，找到最大的Ek
            for k in valid_ecache_list: #迭代所有有效的缓存，找到误差最大的E
                # 不计算k==i节省时间
                if k == i: #不选择和i相等的值
                    continue
                # 计算Ek
                Ek = cal_Ek(os, k)
                # 计算|Ei - Ek|
                delta_E = abs(Ei - Ek)
                # 找到maxDeltaE
                if(delta_E > max_delta_E):
                    max_K = k
                    max_delta_E = delta_E
                    Ej = Ek
            # 返回max_K，Ej
            return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
        # 没有不为0的误差
        else: #第一次循环时是没有有效的缓存值的，所以随机选一个（仅会执行一次）
            # 随机选择alpha_j的索引值
            j = select_j_random(i, os.m)
            # 计算Ej
            Ej = cal_Ek(os, j)
        # 返回j，Ej
        return j, Ej
    
    """
    函数说明：计算Ek,并更新误差缓存
    Parameters:
        os - 数据结构
        k - 标号为k的数据的索引值
    Returns:
        None
    """
    def update_Ek(os, k):
        # 计算Ek
        Ek = cal_Ek(os, k)
        # 更新误差缓存
        os.ecache[k] = [1, Ek]
    
    """
    函数说明：优化的SMO算法
    Parameters:
        i - 标号为i的数据的索引值
        os - 数据结构
    Returns:
        1 - 有任意一对alpha值发生变化
        0 - 没有任意一对alpha值发生变化或变化太小
    """
    def innerL(i, os):
        # 步骤1：计算误差Ei
        Ei = cal_Ek(os, i)
        # 优化alpha,设定一定的容错率
        if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
            # 使用内循环启发方式2选择alpha_j,并计算Ej
            j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
            # 保存更新前的alpha值，使用深层拷贝
            alpha_i_old = os.alphas[i].copy()
            alpha_j_old = os.alphas[j].copy()
            # 步骤2：计算上界H和下界L
            if(os.label_mat[i] != os.label_mat[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.label_mat[j] * (Ei - Ej) / eta
            # 步骤5：修剪alpha_j
            os.alphas[j] = clip_alpha(os.alphas[j], H, L)
            # 更新Ej至误差缓存
            update_Ek(os, j)
            if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
                print("alpha_j变化太小")
                return 0
            # 步骤6：更新alpha_i
            os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
            # 更新Ei至误差缓存
            update_Ek(os, i)
            # 步骤7：更新b_1和b_2:
            b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
            b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
            # 步骤8：根据b_1和b_2更新b
            if(0 < os.alphas[i] < os.C):
                os.b = b1
            elif(0 < os.alphas[j] < os.C):
                os.b = b2
            else:
                os.b = (b1 + b2) / 2.0
            return 1 #表示有更新
        else:
            return 0 #表示没有更新
    
    """
    函数说明：完整的线性SMO算法
    Parameters:
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
        maxIter - 最大迭代次数
    Returns:
        oS.b - SMO算法计算的b
        oS.alphas - SMO算法计算的alphas
    """
    def smo_p(data_mat_in, class_labels, C, toler, max_iter):
        # 初始化数据结构
        os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
        # 初始化当前迭代次数
        iter = 0
        entrie_set = True #是否在全部数据集上迭代
        alpha_pairs_changed = 0
        # 遍历整个数据集alpha都没有更新或者超过最大迭代次数，则退出循环
        while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
            alpha_pairs_changed = 0
            if entrie_set: # 遍历整个数据集
                for i in range(os.m):
                    # 使用优化的SMO算法
                    alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
                    print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
                iter += 1
            # 遍历非边界值
            else:
                # 遍历不在边界0和C的alpha
                non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
                for i in non_bound_i_s:
                    alpha_pairs_changed += innerL(i, os)
                    print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
                iter += 1
            # 遍历一次后改为非边界遍历
            if entrie_set:
                entrie_set = False #进行切换，遍历非边界数据集
            # 如果alpha没有更新，计算全样本遍历
            elif(alpha_pairs_changed == 0):
                entrie_set = True
            print("迭代次数:%d" % iter)
        # 返回SMO算法计算的b和alphas
        return os.b, os.alphas
    
    """
    函数说明：分类结果可视化
    Returns:
        dataMat - 数据矩阵
        classLabels - 数据标签
        w - 直线法向量
        b - 直线截距 
    Returns:
        None
    """
    def show_classifer(data_mat, class_labels, w, b):
        data_mat, label_mat = load_dataset('testSet.txt')
        fig = plt.figure()
        ax = fig.add_subplot(111)
        cm = mpl.colors.ListedColormap(['g', 'b'])
        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()
    """
    函数说明：计算w
    Returns:
        dataArr - 数据矩阵
        classLabels - 数据标签
        alphas - alphas值
    Returns:
        w - 直线法向量
    """
    def cal_w_s(alphas, data_array, class_labels):
        X = np.mat(data_array)
        label_mat = np.mat(class_labels).transpose()
        m, n = np.shape(X)
        w = np.zeros((n, 1))
        for i in range(m):
            w += np.multiply(alphas[i] * label_mat[i], X[i, :].T)
        return w
    
    
    if __name__ == '__main__':
        data_array, class_labels = load_dataset('testSet.txt')
        b, alphas = smo_p(data_array, class_labels, 0.6, 0.001, 40)
        w = cal_w_s(alphas, data_array, class_labels)
        print('b=', b)
        print('w=', w)
        print('alphas=', alphas[alphas > 0])
        for i in range(100):
            if alphas[i] > 0.0:
                print(data_array[i], class_labels[i])
        show_classifer(data_array, class_labels, w, b)

[f03dd12888eb46f6bf1d4dd60dd0dfb5.png]: https://img-blog.csdnimg.cn/f03dd12888eb46f6bf1d4dd60dd0dfb5.png