SMO算法

代码：

import random

import numpy as np

"""
函数说明:随机选择alpha
Parameters:
    i - alpha
    m - alpha参数个数
Returns:
    j -
"""


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


""""
函数说明:根据上下限修剪alpha
Parameters:
    alpha_j - alpha值
    H - alpha上限
    L - alpha下限
Returns:
    alpha_j - alpha值
"""


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


"""
函数说明:简化版SMO算法
Parameters:
    dataMatIn - 数据矩阵 x
    classLabels - 数据标签 y
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
    无
"""

"""
未修剪前
alpha_2_new =alpha_2_old + y_2*(E1-E2)/eta
eta = K_11 + K_22 - 2K_12
alpha_1_new =(C - alpha2_new*y_2)*y_1
"""


def smoSimple(dataMatIn, classLables, C, toler, maxIter):
    dataMatrix = np.mat(dataMatIn)  # 将传入的数据矩阵转换为numpy的mat存储 x
    labelMat = np.mat(classLables).transpose()  # 将传入的数据标签y转置
    b = 0
    m, n = np.shape(dataMatrix)  # 获得数据维度 m行n列，也可以看做有m个数据样本，n个特征
    alphas = np.mat(np.zeros((m, 1)))
    iter_num = 0  # 初始化迭代次数
    # 最多迭代maxIter次
    while (iter_num < maxIter):
        alphaPairsChanged = 0  # 用于记录alphas是否已经进行优化
        for i in range(m):
            # 有m个数据样本，遍历m次
            #  步骤1：计算误差Ei
            # 决策公式f(x) = sign(wx+b) = f(x) = ∑alpha_i*y_i*<x_i,x> + b
            # dataMatrix[i,:]表示矩阵中第i行所有元素
            # np.multiply(a,b) 将a,b中对应元素相乘
            fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            # 求Ei = fXi -y_i
            Ei = fXi - float(labelMat[i])
            # 优化alpha,更设定一定的容错率
            # 下面就是界定非边界的alpha，正负间隔都会被测试，筛选满足KKT条件的，在边界上的alpha等于0或者C，则不再优化
            if ((labelMat[i] * Ei < - toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                # y_i *Ei <toler且alphas[i] <c 或者 y_i*Ei>toler 且alphas[i] >0
                # 随机选择另一个与alphas_i成对优化的alphas_j
                j = selectJrand(i, m)
                # 计算误差Ej
                fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                # 保存更新前的alphas值 在更新或比较时会用到
                # 即alpha1_old ,alpha2_old
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # 步骤2：计算上下界L,H 将数学分析得到的结果直接带入
                if labelMat[i] != labelMat[j]:
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C - alphas[i] + alphas[j])
                else:
                    L = max(0, alphas[i] + alphas[j] - C)
                    H = min(C, alphas[i] + alphas[j])
                if L == H:
                    print('L==H')
                    continue  # 如果上下界相等，则跳出本次for循环

                # 步骤3 计算eta 参数更新里的:K_11 + K_22 - 2K_12,有点不同，是因为假设条件不一样，但是原理相同
                eta = dataMatrix[i, :] * dataMatrix[i, :].T + dataMatrix[j, :] * dataMatrix[j, :].T - 2.0 * dataMatrix[
                                                                                 i, :] * dataMatrix[j, :].T
                if eta <= 0:
                    print('eta <=0')
                    continue  # 满足KKT条件，跳出该循序
                # 步骤4： 否则更新alphas[j]
                alphas[j] = alphas[j] + labelMat[j] * (Ei - Ej) / eta
                # 步骤5：修剪alphas[j]
                alphas[j] = clipAlpha(alphas[j], H, L)
                # abs 函数  返回数字的绝对值
                # 查看更新后的alpha是否满足条件，即变化大不大，如果不大，说明该alpha已经趋于稳定了，不用继续更新了
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    # print(f'对于样本{i},alpha_j变化太小')
                    continue  # 对于第i个样本 alpha_j 变化太小，舍弃这个alpha_i 跳出循环，重新选择第一个变量
                # 若变化满足要求，则更新alpha_i
                # alphas[i] = (C - alphas[j]*labelMat[j])*labelMat[i] =alphas[i] = alphas[i] +(alphaJold - alphas[j])*labelMat[i]*labelMat[j]
                # alpha[i] = alpha[i] + y_j *y_i*(alphaJold - alphas[j])
                alphas[i] = alphas[i] + (alphaJold - alphas[j]) * labelMat[i] * labelMat[j]
                # 更新b
                b1 = b - Ei + (alphaIold - alphas[i]) * labelMat[i] * dataMatrix[i, :] * dataMatrix[i, :].T + (
                        alphaJold - alphas[j]) * labelMat[j] * dataMatrix[j, :] * dataMatrix[i, :].T
                b2 = b - Ej + (alphaIold - alphas[i]) * labelMat[i] * dataMatrix[j, :] * dataMatrix[i, :].T + (
                        alphaJold - alphas[j]) * labelMat[j] * 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  # 当运行到这里说明alpha已经进行更新了，因此更新优化次数加一
                print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            # 若经过了for循环alpha更新次数依然为0
            # 意味着这一轮迭代已经无法再优化参数，故迭代次数加一进入下一轮迭代
            # 若alpha优化次数为0 则进行下一次迭代
            # 若迭代次数一直加到最大迭代次数，则表明次数参数已经接近最优，无法再被优化，则跳出循环
            iter_num += 1
        else:
            # 若alpha更新次数有增加，即alpha还能再优化，则将迭代次数设为0，
            # 将已经更新优化好的参数重新作为起始参数，进入下一轮迭代更新
            iter_num = 0
        print("迭代次数: %d" % iter_num)
        # print('已经跑完一遍')
    return b, alphas


def loadDataSet(filename):
    dataMat = []
    labelMat = []
    fr = open(filename)
    for line in fr.readlines():  # 逐行读取，略除空格
        # print(line)
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])  # 添加数据
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


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