当前位置: article > 正文

【机器学习】支持向量机（6）——SMO算法Python代码实现_svm smo 算法 python 代码机器学习方法李航

作者：Monodyee | 2024-05-06 19:56:13

踩

svm smo 算法 python 代码机器学习方法李航

前言

此文介绍了 $S M O$ 算法，以及前面我们介绍了支持向量机的理论，下面我们就该通过代码来实现了。

由于 $S M O$ 算法不易于理解，为了让大家正确理解它的工作流程，我们先从简化版的 $S M O$ 算法开始讨论。

【数据集以及代码】链接：https://pan.baidu.com/s/1mGQNPA4xmgUEANo1o8rHxg 密码：oq2g

应用简化版 $S M O$ 算法处理小规模数据集

通过之前的学习，我们知道 $S M O$ 算法中的外循环确定要优化的 $\alpha$ ，而简化版的会跳过这一部分，首先在训练数据集上遍历每一个 $\alpha$ ，然后在剩下的 $\alpha$ 集合中随机选取另一个 $\alpha$ ，从而构建 $\alpha$ 对。

这里有一点非常重要，之前我们也提到过，这里再说明一下：
我们需要同时改变两个变量 $\alpha$ ，因为约束条件：
$\alpha_1y_1+\alpha_2y_2 = -\sum_{i=3}^N\alpha_iy_i=\varsigma$
其中， $\varsigma$ 为常数

为此，我们先构建一个辅助函数，用于在某个区间范围内随机选择一个整数；同时我们也需要另一个辅助函数，用于在数值太大时，对其进行调整；还有一个函数是加载数据集的。

代码如下：

import random
from numpy import *

# 加载数据集
def loadDataSet(filename):
    dataMat = []
    labelMat = []
    fr = open(filename)
    for line in fr.readlines(): # 解析文本文件中的数据
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])]) # 数据矩阵
        labelMat.append(float(lineArr[2]))  # 类标签
    return dataMat, labelMat

# 用于在某个区间范围内随机选择一个整数
def selectJrand(i, m):  # i为第一个alpha下标，m是所有alpha对应的数目
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

# 用于调整大于H或小于L的alpha值
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

# 测试
dataArr, labelArr = loadDataSet('testSet.txt')
print(labelArr)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

运行结果1：

[-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0,...
1

可以看出，这里采用的类别标签是-1和+1

简化版 $S M O$ 算法

伪代码如下：

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

python3 代码：

# 简化版SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):  # 数据集、类别标签、常数C、容错率和退出前最大的循环次数
    dataMatrix = mat(dataMatIn)  # 将数据集用mat函数转换为矩阵
    labelMat = mat(classLabels).transpose()  # 转置类别标签（使得类别标签向量的每行元素都和数据矩阵中的每一行一一对应）
    b = 0
    m, n = shape(dataMatrix)  # 通过矩阵的shape属性得到常数m,n
    alphas = mat(zeros((m, 1)))  # 构建一个alpha列矩阵，初始化为0
    iter = 0  # 建立一个iter变量（该变量存储在没有任何alpha改变的情况下遍历数据集的次数，当该变量达到输入值maxIter时，函数结束）
    while (iter < maxIter):  # 遍历数据集的次数小于输入的maxIter
        alphaPairsChanged = 0  # 用于记录alpha是否已经进行优化，先设为0
        for i in range(m):  # 顺序遍历整个集合
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b  # fXi可以计算出来，是我们预测的类别
            Ei = fXi - float(labelMat[i])  # 计算误差，基于这个实例的预测结果和真实结果的比对
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or (
                    (labelMat[i] * Ei > toler) and (alphas[i] > 0)):  # 如果误差过大，那么对该数据实例对应的alpha值进行优化
                j = selectJrand(i, m)  # 随机选取第二个alpha，调用辅助函数selectJrand
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])  # 计算误差（同第一个alpha值的误差计算方法）
                alphaIold = alphas[i].copy()  # 浅复制，两对象互不影响，为了稍后对新旧alpha值进行比较
                alphaJold = alphas[j].copy()
                if (labelMat[i] != labelMat[j]):  # 保证alpha值在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  # 本次循环结束，进入下一次for循环
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - dataMatrix[j, :] * dataMatrix[j, :].T  # alpha值的最优修改值
                if eta >= 0:
                    print("eta >= 0")
                    continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta  # 计算出新的alpha[j]值
                alphas[j] = clipAlpha(alphas[j], H, L)  # 调用clipAlpha辅助函数以及L,H值对其进行调整
                if (abs(alphas[j] - alphaJold) < 0.00001):  # 检查alpha[j]是否有轻微改变，如果是的话，就退出for循环
                    print("J not moving enough!")
                    continue
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])  # 对alpha[i]进行同样的改变，改变大小一样方向相反

                # 给alpha[i]和alpha[j]设置一个常数项b，等同于上篇博客中计算阀值b部分，可以对照看
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[i, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].T
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].T
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2) / 2.0
                alphaPairsChanged += 1  # 当for循环运行到这一行，说明已经成功的改变了一对alpha值，同时让alphaPairsChanged加1
                print("iter:%d i:%d,pairs changed %d" % (iter, i, alphaPairsChanged))
        # for循环结束后，检查alpha值是否做了更新
        if (alphaPairsChanged == 0):  # 如果没有，iter加1
            iter += 1  # 下面后回到while判断
        else:  # 如果有更新则iter设为0后继续运行程序
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas  # （只有在所有数据集上遍历maxIter次，且不再发生任何alpha值修改之后，程序才会停止，并退出while循环）

# 测试
b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
print(b)
print(alphas[alphas > 0])  # 观察元素大于0的

shape(alphas[alphas > 0])  # 得到支持向量的个数

for i in range(100):
    if alphas[i] > 0.0:
        print("简单版支持向量为:", dataArr[i], labelArr[i])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

运行结果2：
这里写图片描述

为了画出图形来，我们还需要计算w的值，如下：

python3 代码：

# 计算w的值
# 该程序最重要的是for循环，for循环中实现的仅仅是多个数的乘积，前面我们计算出的alpha值，大部分是为0
# 虽然遍历数据集中的所有数据，但是起作用的只有支持向量，其他对计算w毫无作用
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


ws = calcWs(alphas, dataArr, labelArr)
print(ws)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

运行结果3：
这里写图片描述

画图代码：

将上面程序运行的结果输入下面代码：

from numpy import *
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

xcord0 = []
ycord0 = []
xcord1 = []
ycord1 = []
markers = []
colors = []
fr = open('testSet.txt')
for line in fr.readlines():
    lineSplit = line.strip().split('\t')
    xPt = float(lineSplit[0])
    yPt = float(lineSplit[1])
    label = int(lineSplit[2])
    if (label == -1):
        xcord0.append(xPt)
        ycord0.append(yPt)
    else:
        xcord1.append(xPt)
        ycord1.append(yPt)

fr.close()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord0, ycord0, marker='s', s=90)
ax.scatter(xcord1, ycord1, marker='o', s=50, c='red')
plt.title('Support Vectors Circled')
circle = Circle((4.658191, 3.507396), 0.5, facecolor='none', edgecolor=(0, 0.8, 0.8), linewidth=3, alpha=0.5)
ax.add_patch(circle)
circle = Circle((3.457096, -0.082216), 0.5, facecolor='none', edgecolor=(0, 0.8, 0.8),linewidth=3, alpha=0.5)
ax.add_patch(circle)
circle = Circle((5.286862, -2.358286), 0.5, facecolor='none', edgecolor=(0, 0.8, 0.8), linewidth=3,alpha=0.5)
ax.add_patch(circle)
circle = Circle((6.080573, 0.418886), 0.5, facecolor='none', edgecolor=(0, 0.8, 0.8), linewidth=3,alpha=0.5)
ax.add_patch(circle)
# plt.plot([2.3,8.5], [-6,6]) #seperating hyperplane
b = -3.82407793
w0 = 0.81085367
w1 = -0.25395222
x = arange(-2.0, 12.0, 0.1)
y = (-w0 * x - b) / w1
ax.plot(x, y)
ax.axis([-2, 12, -8, 6])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

运行结果4：
这里写图片描述

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

下面我们讨论完整版的 $S M O$ 算法。

在这两个版本中，实现 $\alpha$ 的更改和代数运算的优化环节一模一样。在优化过程中，唯一不同的就是选择 $\alpha$ 的方式。完整版 $S M O$ 算法应用了一些能够提速的启发式方法。

$S M O$ 算法是通过一个外循环来选择第一个 $\alpha$ 值的，并且其选择过程会在两种方式之间进行交替：

在所有数据集上进行单遍扫描
在非边界 $\alpha$ 中实现单遍扫描

所谓非边界 $\alpha$ ，即：指那些不等于边界0或C的 $\alpha$ 值。

对整个数据集的扫描相当容易，而实现非边界扫描，首先需要建立这些 $\alpha$ 值的列表，然后再对这个列表进行遍历，同时该步骤会跳过那些已知的不会改变的 $\alpha$ 值。
在选择第一个 $\alpha$ 值后，算法会通过一个内循环来选择第二个 $\alpha$ 值。在优化过程中通过最大化步长来获取第二个 $\alpha$ 值。
我们建立一个全局的缓存用于保存误差值，从中选择使得步长或者 $E_i - E_j$ 最大的 $\alpha$ 值。

利用完整版进行优化

以下为用于清理代码的数据结构和3个用于对 $E$ 进行缓存的辅助函数：

完整版 $S M O$ 的支持函数

python3 代码：

# 完整版Platt SMO的支持函数
# （不加核函数）构建一个仅包含init方法的optStruct类，实现其成员变量的填充
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))

# 计算E值并返回
# 不加核函数
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

# 内循环中的启发式方法，选择第二个alpha
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]  # 构建出一个非零表，函数nonzero返回一个列表

    # 在所有值上进行循环并选择其中使得该变量最大的那个值
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


# 计算误差值并存入缓存中，在对alpha值进行优化之后会用到这个值
def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

完整版$SMO#算法内外循环函数

# 完整Platt SMO算法中的优化例程（与smoSimple()函数一样，只是使用了自己的数据结构，在参数oS中传递）
# innerL()函数用来选择第二个alpha，并在可能是对其进行优化处理，如果存在任意一对alpha值发生变化，那么返回1
# 不加核函数
def innerL(i, oS):
    Ei = calcEk(oS, i)  # 计算E值
    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)):
        j, Ej = selectJ(i, oS, Ei)  # 第二个alpha选择中的启发式方法
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        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
        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
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)  # 计算误差值，并存入缓存
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        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[j, :] * oS.X[j, :].T - oS.labelMat[j] * \
             (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        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

# 完整版Platt SMO的外循环代码
# 不加核函数
def smoP(dataMatIn, classLabels, C, toler, maxIter):  # 数据集、类别标签、常数C、容错率和退出前最大的循环次数
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)  # 构建一个数据结构来容纳所有的数据

    # 初始化一些控制变量
    iter = 0
    entireSet = True
    alphaPairsChanged = 0

    # 代码主体
    # 退出循环的条件：
    # 1、迭代次数超过指定的最大值；
    # 2、历整个集合都没有对任意alpha值进行修改（即：alphaPairsChanged=0）
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:  # 遍历所有的值
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)  # 调用innerL()函数
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
        else:  # 遍历非边界值（非边界值指的是那些不等于边界0或C的alpha值）
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print("non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

测试代码：

# 测试
dataArr, labelArr = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
print(b)
print(alphas[alphas > 0])  # 观察元素大于0的

shape(alphas[alphas > 0])  # 得到支持向量的个数

for i in range(100):
    if alphas[i] > 0.0:
        print("完整版支持向量为:", dataArr[i], labelArr[i])


# 计算w的值
# 该程序最重要的是for循环，for循环中实现的仅仅是多个数的乘积，前面我们计算出的alpha值，大部分是为0
# 虽然遍历数据集中的所有数据，但是起作用的只有支持向量，其他对计算w毫无作用
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


ws = calcWs(alphas, dataArr, labelArr)
print(ws)

# 对数据进行分类处理
dataMat = mat(dataArr)
print(dataMat[0] * mat(ws) + b)
print(labelArr[0])

dataMat = mat(dataArr)
print(dataMat[1] * mat(ws) + b)
print(labelArr[1])

dataMat = mat(dataArr)
print(dataMat[2] * mat(ws) + b)
print(labelArr[2])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

运行结果5：

这里写图片描述

将程序运行结果代入画图代码同样可以画出图：

这里写图片描述

图中我们可以看出，这两个类的数据点分布在一条直线的两侧，我们可以得到两类的分割线。但倘若两类数据集分布在一个圆的内部和外部呢？
之前博客理论部分我们以及说明了非线性支持向量机，使用核技巧

下面我们就来讨论在复杂数据上应用核函数

在复杂数据上应用核函数

代码如下：

import random
from numpy import *

def loadDataSet(filename):
    dataMat = []
    labelMat = []
    fr = open(filename)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# 加核函数
def calcEkK(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek


# 内循环中的启发式方法，选择第二个alpha
def selectJK(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]  # 构建出一个非零表，函数nonzero返回一个列表

    # 在所有值上进行循环并选择其中使得该变量最大的那个值
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEkK(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEkK(oS, j)
    return j, Ej


# 计算误差值并存入缓存中，在对alpha值进行优化之后会用到这个值
def updateEkK(oS, k):
    Ek = calcEkK(oS, k)
    oS.eCache[k] = [1, Ek]


# 核转换函数
def kernelTrans(X, A, kTup):  # 2个数值型变量和一个元祖，kTup是一个包含核函数信息的元祖，该元祖的第一个参数描述的是所用核函数的类型的一个字符串
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    if kTup[0] == 'lin':  # 线性核函数，内积计算在“所有数据集”和“数据集中的一行”这两个输入之间展开
        K = X * A.T
    elif kTup[0] == 'rbf':  # 径向核函数，在for循环中对矩阵的每个元素计算高斯函数值
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        K = exp(K / (-1 * kTup[1] ** 2))
    else:
        raise NameError('Houston We have a  Problem -- That Kernel is not recognized')  # 遇到无法识别的核函数抛出异常
    return K


class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)


# 加入核函数
def innerLK(i, oS):
    Ei = calcEkK(oS, i)  # 计算E值
    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)):
        j, Ej = selectJK(i, oS, Ei)  # 第二个alpha选择中的启发式方法
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        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
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
        if eta >= 0:
            print("eta>=0")
            return 0
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEkK(oS, j)  # 计算误差值，并存入缓存
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        updateEkK(oS, i)
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * \
             (oS.alphas[j] - alphaJold) * oS.K[i, j]
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * \
             (oS.alphas[j] - alphaJold) * oS.K[j, j]
        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 smoPK(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerLK(i, oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerLK(i, oS)
                print("non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False  # toggle entire set loop
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas


# 利用核函数进行分类的径向基测试函数
def testRbf(k1=1.3):  # 参数为高斯径向基函数中的一个用户定义变量
    dataArr, labelArr = loadDataSet('testSetRBF.txt')  # 读入数据集
    b, alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))  # 调用函数smoP,核函数类型为"rbf"
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]  # 构建支持向量矩阵
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0

    for i in range(m):
        # 整个代码最重要部分（for循环开始的这两行），给出如何利用核函数进行分类
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))  # 调用kernelTrans()核转换函数
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b  # 与其面的alpha及类别标签值求积
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):  # 与前面一个for循环相比，只是数据集不同，使用的是测试数据集
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))


# 测试
testRbf()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

运行结果6：
这里写图片描述

你可以尝试更换不同的 $k 1$ 值以观察测试错位率、训练错误率，支持向量个数

参考书籍《机器学习实战》

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Monodyee/article/detail/545711

【机器学习】支持向量机（6）——SMO算法Python代码实现_svm smo 算法 python 代码 机器学习方法 李航

前言

应用简化版 S M O SMO SMO算法处理小规模数据集

利用完整版进行优化

在复杂数据上应用核函数

【机器学习】支持向量机（6）——SMO算法Python代码实现_svm smo 算法 python 代码机器学习方法李航

应用简化版 $S M O$ 算法处理小规模数据集