SMO算法在软间隔SVM分类中的应用

序列最小优化算法（SMO）

SMO算法是一种启发式的算法，它在求解过程中通过以分析的方式来定位最优解可能存在的位置，从而避免了传统方法在求解中所遭遇的大量数值计算问题，并且最终以迭代的方式来求得最优解。

SMO算法原理

获得没有修剪的原始解；

对原始解进行修建；

更新阈值b

#SMO算法在软间隔SVM分类中的应用
import numpy as np
import random
import matplotlib.pyplot as plt

def simple_smo(dataset, labels, C, max_iter):
    ''' 简化版SMO算法实现，未使用启发式方法对alpha对进行选择.
    :param dataset: 所有特征数据向量
    :param labels: 所有的数据标签
    :param C: 软间隔常数, 0 <= alpha_i <= C
    :param max_iter: 外层循环最大迭代次数
    '''
    dataset = np.array(dataset)
    m, n = dataset.shape #样本数量，特征数量
    labels = np.array(labels)

    # 初始化参数λ，b为0
    lambds = np.zeros(m) #每个样本都有一个λ乘子
    b = 0
    it = 0


    while it < max_iter:
        pair_changed = 0 #选取的一对值相较于之前是否有变化，没有的话可以停止迭代
        for i in range(m):
            λ_i, x_i, y_i = lambds[i], dataset[i], labels[i] #选取一组λ
            fx_i = SVM_predict(x_i,lambds,dataset,labels,b)
            E_i = fx_i - y_i
            j = select_j(i, m) #选取另一个λ
            λ_j, x_j, y_j = lambds[j], dataset[j], labels[j]
            fx_j = SVM_predict(x_j,lambds,dataset,labels,b)
            E_j = fx_j - y_j
            K_ii, K_jj, K_ij = np.dot(x_i, x_i), np.dot(x_j, x_j), np.dot(x_i, x_j)
            eta = K_ii + K_jj - 2*K_ij #
            if eta <= 0:
                print('WARNING  eta <= 0')
                continue
            # 获取更新的alpha对
            λ_i_old, λ_j_old = λ_i, λ_j #未更新前的参数
            λ_j_new = λ_j_old + y_j*(E_i - E_j)/eta
            # 对alpha进行修剪
            if y_i != y_j:
                L = max(0, λ_j_old - λ_i_old)
                H = min(C, C + λ_j_old - λ_i_old)
            else:
                L = max(0, λ_i_old + λ_j_old - C)
                H = min(C, λ_j_old + λ_i_old)
            λ_j_new = clip(λ_j_new, L, H) #根据上下界修剪
            λ_i_new = λ_i_old + y_i*y_j*(λ_j_old - λ_j_new) #根据公式反推另一个参数
            if abs(λ_j_new - λ_j_old) < 0.00001: #这个参数已经优化到最佳，换下一个
                #print('WARNING   alpha_j not moving enough')
                continue

            #更新b
            lambds[i], lambds[j] = λ_i_new, λ_j_new
            b_i = -E_i - y_i*K_ii*(λ_i_new - λ_i_old) - y_j*K_ij*(λ_j_new - λ_j_old) + b
            b_j = -E_j - y_i*K_ij*(λ_i_new - λ_i_old) - y_j*K_jj*(λ_j_new - λ_j_old) + b
            if 0 < λ_i_new < C:
                b = b_i
            elif 0 < λ_j_new < C:
                b = b_j
            else:
                b = (b_i + b_j)/2
            pair_changed += 1
            print('INFO   iteration:{}  i:{}  pair_changed:{}'.format(it, i, pair_changed))
        if pair_changed == 0: #参数优化完成，下一轮
            it += 1
        else: #参数没有优化完成，继续迭代
            it = 0
        print('iteration number: {}'.format(it))
    return lambds, b

def SVM_predict(x,lambds,data,label,b):
    "SVM分类器函数 y = w^Tx + b,即文中的f(x)"
    res = 0
    for i in range(data.shape[0]):
        res += lambds[i]*label[i]*(data[i].dot(x.T))
    return res + b

def get_w(lambdas, dataset, labels):
    #通过λ求w
    w = 0
    for i in range(len(dataset)):
        w += lambdas[i]*y[i]*dataset[i]
    return w

def clip(alpha, L, H):
    #修建λ的值到L和H之间.

    if alpha < L:
        return L
    elif alpha > H:
        return H
    else:
        return alpha

def select_j(i, m):
    #在m中随机选择除了i之外剩余的
    l = list(range(m))
    seq = l[: i] + l[i+1:]
    return random.choice(seq)

def get_point():
    x_true =  [[1,1.5],[1,3],[4,5],[2,4]]
    x_false = [[1,0.5],[4,2],[5,1],[4,1]]
    x_true =  [[1,2],[1,3],[4,5],[2,4],[2,1.5]]
    x_false = [[1,0.5],[4,2],[5,1],[4,1],[3,2]]
    x_all = np.array(x_true+x_false)
    y = [1]*len(x_true) + [-1]*len(x_false)
    return x_all,y,x_true,x_false

def plot(x_true,x_false,w,b):
    plot_x = np.arange(0,7,0.1)
    plot_y = -(w[0]*plot_x+b)/w[1]
    plt.scatter([x[0] for x in x_true],[x[1] for x in x_true] , c='r' , label='+1')
    plt.scatter([x[0] for x in x_false],[x[1] for x in x_false] , c='b',label='-1')
    plt.plot(plot_x,plot_y,c = 'green')
    plt.xlim(0,7)
    plt.ylim(0,7)
    plt.legend()
    plt.plot()
    plt.show()

if __name__ == '__main__':
    x,y,x_true,x_false = get_point()
    MAX_EPOCHES = 20 #最大迭代次数
    C = 50 #λ边界
    #调节C的值会有不同的改变
    lambdas, b = simple_smo(x, y, C, MAX_EPOCHES)
    w = get_w(lambdas, x, y)
    print('-'*40+'result'+'-'*40)
    print('lambdas:{}\nw:{}\nb:{}'.format(lambdas,w,b))
    plot(x_true,x_false,w,b) #LHQ