绘制出简化的SMO算法和完整SMO算法在testSet.txt数据集上的分类超平面

首先介绍一下简化版SMO算法（Sequential Minimal Optimization），其主要思想是每次选择两个变量进行优化，而其他变量保持不变。简化版SMO算法的步骤如下：

1. 初始化参数 alpha、b 和误差缓存 E
2. 选择两个不同的变量 alpha1 和 alpha2
3. 计算这两个变量对应样本的误差 E1 和 E2
4. 计算变量 alpha2 的上下界 L 和 H
5. 如果 L 和 H 相等，不做任何处理直接退出本轮循环
6. 计算变量 alpha2 的新值
7. 如果新值变化太小，不做任何处理直接退出本轮循环
8. 计算变量 alpha1 的新值
9. 更新参数 b 和误差缓存 E
10. 如果有任意 alpha 值发生变化，增加已经遍历的 alpha 对数目，否则继续循环

下面是完整版SMO算法的步骤：

1. 初始化参数 alpha、b 和误差缓存 E
2. 选择两个不同的变量 alpha1 和 alpha2
3. 计算这两个变量对应样本的误差 E1 和 E2
4. 计算变量 alpha2 的上下界 L 和 H
5. 如果 L 和 H 相等，不做任何处理直接退出本轮循环
6. 计算变量 alpha2 的新值
7. 如果新值变化太小，不做任何处理直接退出本轮循环
8. 计算变量 alpha1 的新值
9. 更新参数 b 和误差缓存 E
10. 如果有任意 alpha 值发生变化，增加已经遍历的 alpha 对数目，否则继续循环
11. 在所有非边界 alpha 中选择违反 KKT 条件的变量作为 alpha1
12. 如果没有这样的变量，随机选择一个 alpha1
13. 选择 alpha2，使得优化步骤最大
14. 进行步骤 3-10，直到达到最大迭代次数或者所有 alpha 值都不再变化

下面是在testSet.txt数据集上绘制分类超平面的代码实现：

import numpy as np
import matplotlib.pyplot as plt

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(np.random.uniform(0, m))
    return j

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

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    b = 0
    m, n = np.shape(dataMatrix)
    alphas = np.mat(np.zeros((m, 1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i, m)
                fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                if (labelMat[i] != labelMat[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
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T
                if eta >= 0:
                    print("eta >= 0")
                    continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    print("j not moving enough")
                    continue
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                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
                print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas

def plotSVM(dataMatIn, classLabels, b, alphas):
    dataArr = np.array(dataMatIn)
    n = np.shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(classLabels[i]) == 1:
            xcord1.append(dataArr[i, 0])
            ycord1.append(dataArr[i, 1])
        else:
            xcord2.append(dataArr[i, 0])
            ycord2.append(dataArr[i, 1])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    for i in range(n):
        if alphas[i] > 0:
            ax.plot(dataArr[i, 0], dataArr[i, 1], 'o')
    ws = np.zeros((2, 1))
    for i in range(n):
        ws += np.multiply(alphas[i] * classLabels[i], dataMatrix[i, :].T)
    x = np.arange(-2.0, 12.0, 0.1)
    y = (-b - ws[0] * x) / ws[1]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet('testSet.txt')
    b, alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    plotSVM(dataMat, labelMat, b, alphas)

运行上述代码后，即可在图形界面上看到绘制的分类超平面。