机器学习之支持向量机(SVM)的求解方法

Posted Icy Hunter

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了机器学习之支持向量机(SVM)的求解方法相关的知识,希望对你有一定的参考价值。

文章目录

前言

支持向量机就是寻找一个超平面,将不同的样本分分隔开来,其中间隔分为硬间隔和软间隔,硬间隔就是不允许样本分错,而软间隔就是允许一定程度上样本存在偏差,后者更符合实际。

支持向量机思路简单但是求解过程还是比较复杂,需要将原函数通过拉格朗日乘子法并附上KKT条件是的问题有强对偶性,再使用SMO等算法进行高效的求解。
推导过程可以参考:
机器学习之支持向量机之线性可分型原理介绍及代码实现(SVM)
下面主要实现模型的求解方法。

梯度下降法

梯度下降法是一种比较普适的方法,当模型无法得出解析解,或者解析解求解困难的时候,都可以使用梯度下降法来近似求解。因为梯度下降法需要一轮轮迭代,也需要定义损失函数,因此一般而言,梯度下降法只能获得近似最优解。
SVM可以使用梯度下降法求解,不过得出的解大概率只是近似解,并且不一定满足SVM的公式里的约束条件。

# coding=utf8
import sys;
import random;
import numpy as np
import math

EPS = 0.000000001  # 很小的数字,用于判断浮点数是否等于0
def load_data(filename, dim):
    '''
    输入数据格式: label\\tindex1:value1\\tindex2:value2\\tindex3:value3..., 其中index是特征的编号, 从1开始
    data的数据格式: [[label, sample],[label, sample], ...],  其中sample: [v0, v1, v2, v3, ..., v[dim]]
    '''
    label_ = []
    data_ = []
    for line in open(filename, 'rt'):
        sample = [0.0 for v in range(0, dim + 1)];
        line = line.rstrip("\\r\\n\\t ");
        fields = line.split("\\t");
        label = int(fields[0]);  # LABEL取值: 1 or -1
        sample[0] = 1.0;  # sample第一个元素用于存x0特征, 默认置为1.0[方便把 WX+b => WX]
        for field in fields[1:]:
            kv = field.split(":");
            idx = int(kv[0]);  # ensure idx >= 1
            val = float(kv[1]);
            sample[idx] = val;
        label_.append(label)
        data_.append(sample)
    
    label_ = np.array(label_)
    data_ = np.array(data_)
    return label_, data_
    


def svm_train(train_x, train_y, dim, iterations, lm, lr):
    '''
    data4train:数据集
    dim:样本特征维度
    W:SVM模型的权重
    iterations:迭代次数
    目标函数: obj(<X,y>, W) = (对所有<X,y>SUMmax0, 1 - W*X*y) + lm / 2 * ||W||^2, 即:hinge+L2
    '''
    X = np.zeros(dim + 1)  # <sample, label> => <X, y>
    grad = np.zeros(dim + 1) # 梯度
    num_train = len(train_x);
    global W
    for i in range(0, iterations):
        # 每次迭代随机选择一个训练样本
        index = random.randint(0, num_train - 1);
        y = train_y[index]  # y其实就是label
        X = train_x[index]

        # 计算梯度
        # for j in range(0, dim + 1):
        #	grad = lm * W[j];
        WX = 0.0
        WX += (W * X).sum() 

        if 1 - WX * y > 0:
            grad = lm * W - X * y
        else:  # 1-WX *y <= 0的时候,目标函数的前半部分恒等于0, 梯度也是0
            grad = lm * W - 0;
        # 更新权重, lr是学习速率
        W = W - lr * grad
        


def svm_predict(x, y, dim, W):
    num_test = len(x);
    num_correct = 0;
    for i in range(0, num_test):
        target = y[i]  # 即label
        X = x[i]  # 即sample
        sum = 0.0;
        sum += (X * W).sum()
        predict = -1;
        # print sum;
        if sum > 0:  # 权值>0,认为目标值为1
            predict = 1;
        if predict * target > 0:  # 预测值和目标值符号相同
            num_correct += 1;

    return num_correct * 1.0 / num_test;


if __name__ == "__main__":
    # 设置参数
    epochs = 10;  # 迭代轮数
    iterations = 10;  # 每一轮中梯度下降迭代次数, 这个其实可以和epochs合并为一个参数
    lm = 0.0001;  # lambda, 对权值做正则化限制的权重
    lr = 0.01;  # lr, 是学习速率,用于调整训练收敛的速度
    dim = 1000;  # dim, 特征的最大维度, 所有样本不同特征总数
    W = np.zeros(dim + 1)  # 权值
    # 导入测试集&训练集
    train_y, train_x = load_data("train.txt", dim)
    test_y, test_x = load_data("test.txt", dim)

    # # 训练, 实际迭代次数=epochs * iterations
    for i in range(0, epochs):
        svm_train(train_x, train_y, dim, iterations, lm, lr);
        accuracy = svm_predict(test_x, test_y, dim, W);

        print("epoch:%d\\taccuracy:%f" % (i, accuracy));
    # 输出结果权值
    for i in range(0, dim + 1):
        if math.fabs(W[i]) > EPS:
            print("权值W%d\\t%f" % (i, W[i]));
    print(W)

这个是参考支持向量机SVM-手写笔记&手动实现这篇博客的,代码改成numpy进行运算了,进行改动的过程也是读懂代码的过程。

这个例子使用的梯度下降法,损失函数应该大概也许相当于用了max(0, x)。整体思路就是,每次就选一个样本点进行参数更新,如果这个样本点对于当前的参数能够正确分类,那么就不更新,如果不能正确分类,就更新。

运行结果如下:

因为每次随机选择样本,那么其实找到的那个超平面大概率不会将样本集完全分开的,但是事实证明梯度下降法还是有效果的,所有,应该可以勉强认为训练出来的是软间隔的SVM吧?

SMO算法

贴上b站视频,以便日后再去看,实话说目前也没有完全走通这个算法的流程。
快速理解SMO算法
SMO算法思路很简单,因为存在约束条件

所以每次更新两个α,剩余看成常量,α能够通过看成的常量导出,然后满足约束条件的情况下求出极值,更新α。每次都更新两个,当固定其他的α时,能够求出选取α的更新的解析解,所以就算起来非常快。
那么该如何选择呢?
第一个αi应该选择违反KKT条件最大的。
第二个αj应该选择于第一个αi差值最大的。
这样能够保证每次更新都是向最快的方向进行更新。
思路是很简单,但是实践起来还是很困难的,因为里面涉及到许多的约束条件,不同情况下需要分类讨论等等。
直接贴上别人的代码吧

# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
import random
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 
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 print_m(self):
        print("self.X", self.X)
        print("self.labelMat", self.labelMat)
        print("self.C", self.C)
        print("self.tol", self.tol)
        print("self.m", self.m)
        print("self.alphas", self.alphas)
        print("self.b", self.b)
        print("self.eCache", self.eCache)

# 随机选取一个J值
def selectJrand(i,m): 
    j=i
    while (j==i):
        j=int(random.uniform(0,m))
    return j
 
# 根据关于α_1与α_2的优化问题对应的约束问题分析,对α进行截取约束
# 保证取值范围符合约束条件 
def clipAlpha(aj,H,L):  
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj
        
#计算每个样本点k的Ek值,就是计算误差值=预测值-标签值        
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
 

#内循环的启发式方法,获取最大差值|Ei-Ej|对应的Ej的索引J
def selectJ(i, oS, Ei):
    # 传入第一个α对应的索引i和误差值Ei

    maxK = -1 #用于保存临时最大索引
    maxDeltaE = 0 #用于保存临时最大差值--->|Ei-Ej|
    Ej = 0 #保存我们需要的Ej误差值
    oS.eCache[i] = [1,Ei]
    # 获取有效误差的下标
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]  
    if (len(validEcacheList)) > 1: # 如果存在符合条件的αj
        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 # 找到与αi差值最大的αj(最优更新方向)
    else: # 不存在就任意选取
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

# 更新Ek
def updateEk(oS, k): 
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
 
# 实现内循环函数
def innerL(i, oS):
    # 确定了第一个αi
    # 计算出Ei
    Ei = calcEk(oS, i) 
    # 如果违背了KKT条件
    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)): 
        # 选择i对应确定的αj
        j,Ej = selectJ(i, oS, Ei)
        # 因为后面还需要用到原始值,因此需要记录下来
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        # 确定αj的取值范围
        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
        # 2*xi*xj - xi^2 - xj^2 > 0 <=>
        # xi^2 + xj^2 - 2xi*xj < 0  <=>
        # (xi - xj)^2 < 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
        # 满足条件后更新αj
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta 
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        # 更新Ek 
        updateEk(oS, j)
        # 检测变化量是否显著
        if (abs(oS.alphas[j] - alphaJold) < oS.tol): 
            print("j not moving enough")
            return 0
        # αj 的变化量显著则更新 αi
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        # 更新Ek
        updateEk(oS, i)
        # 更新b 
        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
        #根据统计学习方法中阈值b在每一步中都会进行更新,
        #1.当新值alpha_1不在界上时(0<alpha_1<C),b_new的计算规则为:b_new=b1
        #2.当新值alpha_2不在界上时(0 < alpha_2 < C),b_new的计算规则为:b_new = b2
        #3.否则当alpha_1和alpha_2都不在界上时,b_new = 1/2(b1+b2)
        if (0 < oS.alphas[i]<oS.C):
            oS.b = b1
        elif (0 < oS.alphas[j]<oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        # 更新完毕
        return 1
    else: # 没有违反KKT条件
        return 0
 
 
# 根据西瓜书6.37计算W参数
def calcWs(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = array(alphas), array(dataMat), array(labelMat)
    w = dot((tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()
 
 
def smoP(dataMatIn, classLabels, C, toler, maxIter): 
    # 初始化数据
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    # oS.print_m()
    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 += innerL(i,oS) # innerL成功更新return 1
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) 
            iter += 1
        else:
            # 获取非边界值得索引
            nonBoundIs = nonzero((oS.alphas.A > 以上是关于机器学习之支持向量机(SVM)的求解方法的主要内容,如果未能解决你的问题,请参考以下文章

机器学习之SVM支持向量机

机器学习之支持向量机(手推公式版)

机器学习之支持向量机(手推公式版)

机器学习之支持向量机

python机器学习之SVM(支持向量机)实例

机器学习之SVM(支持向量机)