Educoder 机器学习 SMO进行优化 第1关：SMO高效优化算法

第1关：SMO高效优化算法

• 任务要求
• 参考答案
• 评论1

• 任务描述
• 相关知识
  • 了解SMO的工作原理
• 编程要求
• 测试说明

---

任务描述

本关任务：编写一个SMO算法。

相关知识

为了完成本关任务，你需要掌握：1.了解SMO的工作原理，2.如何编写SMO算法

了解SMO的工作原理

SMO算法的目标是求出一系列α和b,一旦求出了这些α, 就很容易计算出权重向量w并得到分隔超平面。 SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的α,那么就增大其中一个同时减小另一个。这里所谓的“合适” 就是指两个α必须要符合一定的条件，条件之一就是这两个α必须要在间隔边界之外，而其第二个条件则是这两个α还没有进行过区间化处理或者不在边界上。 应用简化版SMO算法处理小规模数据集 简化版代码虽然量少,但执行速度慢。Platt SMO算法中的外循环确定要优化的最佳α对。而简化版却会跳过这一部分，首先在数据集上遍历每一个α , 然后在剩下的α集合中随机选择另一个α,从而构建α对。这里有一点相当重要，就是我们要同时改变两个α 。之所以这样做是因为我们有一个约束条件：

为此，我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们也需要另一个辅助函数，用于在数值太大时对其进行调整。

1. Parameters:
2. fileName - 文件名
3. Returns:
4. dataMat - 数据矩阵
5. labelMat - 数据标签
6. def loadDataSet(fileName):
7. dataMat = []; labelMat = []
8. fr = open(fileName)
9. for line in fr.readlines(): #逐行读取，滤除空格等
10. lineArr = line.strip().split('\t')
11. dataMat.append([float(lineArr[0]), float(lineArr[1])]) #添加数据
12. labelMat.append(float(lineArr[2])) #添加标签
13. return dataMat,labelMat
14. ```python
15. """
16. 函数说明:随机选择alpha
18. Parameters:
19. i - alpha
20. m - alpha参数个数
21. Returns:
22. j -
23. """
24. def selectJrand(i, m):
25. j = i #选择一个不等于i的j
26. while (j == i):
27. j = int(random.uniform(0, m))
28. return j
30. """

第一个函数就是我们所熟知的loadDatSet()函数，该函数打开文件并对其进行逐行解析，从而得到每行的类标签和整个数据矩阵。 下一个函数seelectJrand()有两个参数值，其中i是第一个α的下标，m是所有α的数目。只要函数值不等于输人值i，函数就会进行随机选择。 简化版的SMO函数如下： 创建一个alpha向量并将其初始化为O向量 当迭代次数小于最大迭代次数时（外循环） 对数据集中的每个数据向量（内循环）： 如果该数据向量可以被优化： 随机选择另外一个数据向量 同时优化这两个向量 如果两个向量都不能被优化，退出内循环 如果所有向量都没被优化，增加迭代数目，继续下一次循环 完整代码如下：

1. """
2. 函数说明:简化版SMO算法
4. Parameters:
5. dataMatIn - 数据矩阵
6. classLabels - 数据标签
7. C - 松弛变量
8. toler - 容错率
9. maxIter - 最大迭代次数
10. Returns:
11. 无
12. """
13. def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
14. #转换为numpy的mat存储
15. dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
16. #初始化b参数，统计dataMatrix的维度
17. b = 0; m,n = np.shape(dataMatrix)
18. #初始化alpha参数，设为0
19. alphas = np.mat(np.zeros((m,1)))
20. #初始化迭代次数
21. iter_num = 0
22. #最多迭代matIter次
23. while (iter_num < maxIter):
24. alphaPairsChanged = 0
25. for i in range(m):
26. #步骤1：计算误差Ei
27. fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
28. Ei = fXi - float(labelMat[i])
29. #优化alpha，更设定一定的容错率。
30. if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
31. #随机选择另一个与alpha_i成对优化的alpha_j
32. j = selectJrand(i,m)
33. #步骤1：计算误差Ej
34. fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
35. Ej = fXj - float(labelMat[j])
36. #保存更新前的aplpha值，使用深拷贝
37. alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
38. #步骤2：计算上下界L和H
39. if (labelMat[i] != labelMat[j]):
40. L = max(0, alphas[j] - alphas[i])
41. H = min(C, C + alphas[j] - alphas[i])
42. else:
43. L = max(0, alphas[j] + alphas[i] - C)
44. H = min(C, alphas[j] + alphas[i])
45. if L==H: print("L==H"); continue
46. #步骤3：计算eta
47. eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
48. if eta >= 0: print("eta>=0"); continue
49. #步骤4：更新alpha_j
50. alphas[j] -= labelMat[j]*(Ei - Ej)/eta
51. #步骤5：修剪alpha_j
52. alphas[j] = clipAlpha(alphas[j],H,L)
53. if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
54. #步骤6：更新alpha_i
55. alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
56. #步骤7：更新b_1和b_2
57. b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
58. b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
59. #步骤8：根据b_1和b_2更新b
60. if (0 < alphas[i]) and (C > alphas[i]): b = b1
61. elif (0 < alphas[j]) and (C > alphas[j]): b = b2
62. else: b = (b1 + b2)/2.0
63. #统计优化次数
64. alphaPairsChanged += 1
65. #打印统计信息
66. print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
67. #更新迭代次数
68. if (alphaPairsChanged == 0): iter_num += 1
69. else: iter_num = 0
70. print("迭代次数: %d" % iter_num)
71. return b,alphas
73. """

在有几百个点组成的小规模数据集上，简化版的SMO算法的运行时没有什么问题的，但是在更大的数据集上的运行速度就会变慢。在这两个版本中，实现α的更改和代数运算的优化环节一模一样。在优化过程中，唯一的不同就是选择α的方式。完整版的Platt SMO算法应用了一些能够提速的启发方法。 Platt SMO算法是通过一个外循环来选择第一个α值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界α中实现单遍扫描。而所谓非边界α指的就是那些不等于边界0或α值。对整个数据集的扫描相当容易，而实现非边界α值的扫描时，首先需要建立这些α值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的α值。 在选择第一个α值后，算法会通过一个内循环来选择第二个α值。在优化过程中，会通过最大化步长的方式来获得第二个α值。在简化版SMO算法中，我们会在选择j之后计算错 误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的α值。

编程要求

根据提示，在右侧编辑器补充代码，编写完整版的SMO算法，

测试说明

平台会对你编写的代码进行测试：

开始你的任务吧，祝你成功！

import numpy as np
import random
def calcEk(oS, k):
    """
    计算误差
    Parameters：
        oS - 数据结构
        k - 标号为k的数据
    Returns:
        Ek - 标号为k的数据误差
    """
    fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def loadDataSet(fileName):
    """
    读取数据
    Parameters:
        fileName - 文件名
    Returns:
        dataMat - 数据矩阵
        labelMat - 数据标签
    """
    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 selectJrand(i, m):
    """
    函数说明:随机选择alpha_j的索引值

    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

def selectJ(i, oS, Ei):
    """
    内循环启发方式2
    Parameters：
        i - 标号为i的数据的索引值
        oS - 数据结构
        Ei - 标号为i的数据误差
    Returns:
        j, maxK - 标号为j或maxK的数据的索引值
        Ej - 标号为j的数据误差
    """
    maxK = -1; maxDeltaE = 0; Ej = 0                         #初始化
    oS.eCache[i] = [1,Ei]                                      #根据Ei更新误差缓存
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]        #返回误差不为0的数据的索引值
    if (len(validEcacheList)) > 1:                            #有不为0的误差
        for k in validEcacheList:                           #遍历,找到最大的Ek
            if k == i: continue                             #不计算i,浪费时间
            Ek = calcEk(oS, k)                                #计算Ek
            deltaE = abs(Ei - Ek)                            #计算|Ei-Ek|
            if (deltaE > maxDeltaE):                        #找到maxDeltaE
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej                                        #返回maxK,Ej
    else:                                                   #没有不为0的误差
        j = selectJrand(i, oS.m)                            #随机选择alpha_j的索引值
        Ej = calcEk(oS, j)                                    #计算Ej
    return j, Ej                                             #j,Ej

def updateEk(oS, k):
    """
    计算Ek,并更新误差缓存
    Parameters：
        oS - 数据结构
        k - 标号为k的数据的索引值
    Returns:
        无
    """
    Ek = calcEk(oS, k)                                        #计算Ek
    oS.eCache[k] = [1,Ek]                                    #更新误差缓存


def clipAlpha(aj,H,L):
    """
    修剪alpha_j
    Parameters:
        aj - alpha_j的值
        H - alpha上限
        L - alpha下限
    Returns:
        aj - 修剪后的alpah_j的值
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

class optStruct:
    """
    数据结构，维护所有需要操作的值
    Parameters：
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
    """
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn                                #数据矩阵
        self.labelMat = classLabels                        #数据标签
        self.C = C                                         #松弛变量
        self.tol = toler                                 #容错率
        self.m = np.shape(dataMatIn)[0]                 #数据矩阵行数
        self.alphas = np.mat(np.zeros((self.m,1)))         #根据矩阵行数初始化alpha参数为0
        self.b = 0                                         #初始化b参数为0
        self.eCache = np.mat(np.zeros((self.m,2)))         #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。

def innerL(i, oS):
    """
    优化的SMO算法
    Parameters：
        i - 标号为i的数据的索引值
        oS - 数据结构
    Returns:
        1 - 有任意一对alpha值发生变化
        0 - 没有任意一对alpha值发生变化或变化太小
    """
    #步骤1：计算误差Ei
    Ei = calcEk(oS, i)
    #优化alpha,设定一定的容错率。
    if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
        #使用内循环启发方式2选择alpha_j,并计算Ej
        j,Ej = selectJ(i, oS, Ei)
        #保存更新前的aplpha值，使用深拷贝
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        #步骤2：计算上下界L和H
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L==H")
            return 0
        #步骤3：计算eta
        eta = 2.0 * oS.X[i,:] * oS.X[j,:].T - oS.X[i,:] * oS.X[i,:].T - oS.X[j,:] * oS.X[j,:].T
        if eta >= 0:
            print("eta>=0")
            return 0
        #步骤4：更新alpha_j
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
        #步骤5：修剪alpha_j
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        #更新Ej至误差缓存
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("alpha_j变化太小")
            return 0
        #步骤6：更新alpha_i
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        #更新Ei至误差缓存
        updateEk(oS, i)
        #步骤7：更新b_1和b_2
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        #步骤8：根据b_1和b_2更新b
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0
def smoP(dataMatIn, classLabels, C, toler, maxIter):
    """
    完整的线性SMO算法
    Parameters：
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
        maxIter - 最大迭代次数
    Returns:
        oS.b - SMO算法计算的b
        oS.alphas - SMO算法计算的alphas
    """
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)                    #初始化数据结构
    iter = 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 = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]                        #遍历不在边界0和C的alpha
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:                                                                                #遍历一次后改为非边界遍历
            entireSet = False
        elif (alphaPairsChanged == 0):                                                                #如果alpha没有更新,计算全样本遍历
            entireSet = True
        print("迭代次数: %d" % iter)
    return oS.b,oS.alphas    
                                                                    #返回SMO算法计算的b和alphas
def calcWs(alphas,dataArr,classLabels):
    """
    计算w
    Parameters:
        dataArr - 数据矩阵
        classLabels - 数据标签
        alphas - alphas值
    Returns:
        w - 计算得到的w
    """
    X = np.mat(dataArr); labelMat = np.mat(classLabels).transpose()
    m,n = np.shape(X)
    w = np.zeros((n,1))
    for i in range(m):
        w += np.multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

if __name__ == '__main__':
    dataArr, classLabels = loadDataSet('./src/step2/testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)