支持向量机原理介绍

Posted 2021-04-26 一只写诗的程序员

我们之前已经推导出来支持向量机的数学公式，现在我们打算用SMO算法来进行求解，得到这个公式的最佳系数。

我这里先介绍这个算法的思想，介绍完思想后，我们再将思想套入进来。由于目标函数为凸函数（这里凸函数有兴趣的同学可以参考一下数学分析中的凸函数，凸函数等价于这个函数的导函数是单调递增的函数，深入的研究也很多性质，比如詹森不等式这样的）我们可以通过梯度的方法进行优化，但是目标函数是有条件限制的，所以，为了适应梯度个更新，我们更新时一次更新两个变量，我们先列出算法步骤，再一步步讲解这个算法。这里借用维基百科上一般步骤：

首先我们需要进行数学推导以上算法:

我们得到该推导式就可以用代码实现了（注该代码参考自机器学习实战，并参考自知乎文章https://zhuanlan.zhihu.com/p/29604517 这里由于代码存在可读性比较差，这里略做修改）：

# -*- coding:UTF-8 -*- from time import sleep import matplotlib.pyplot as plt import numpy as np import random import types   def loadDataSet(fileName):     '''     读取数据，转换成数据矩阵和结果标签矩阵     :param fileName: 文件     :return:     '''     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):     '''     随机选择alpha     :param i:alpha     :param m:alpha参数个数     :return:     '''     j = i                                 #选择一个不等于i的j     while (j == i):         j = int(random.uniform(0, m))     return j   def clipAlpha(aj,H,L):     '''     修剪alpha     :param aj:alpha值     :param H:alpha上限     :param L:alpha下限     :return:alpah值     '''     if aj > H:         aj = H     if L > aj:         aj = L     return aj     def smoSimple(dataMatIn, classLabels, C, toler, maxIter):     '''     简化版SMO算法     :param dataMatIn: 数据矩阵     :param classLabels: 数据标签     :param C: 松弛变量     :param toler: 容错率     :param maxIter: 最大迭代次数     :return:     '''     #转换为numpy的mat存储     dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()     #初始化b参数，统计dataMatrix的维度     b = 0; m,n = np.shape(dataMatrix)     #初始化alpha参数，设为0     alphas = np.mat(np.zeros((m,1)))     #初始化迭代次数     iter_num = 0     #最多迭代matIter次     while (iter_num < maxIter):         alphaPairsChanged = 0         for i in range(m):             #步骤1：计算误差Ei             fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b             Ei = fXi - float(labelMat[i])             #优化alpha，更设定一定的容错率。             if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):                 #随机选择另一个与alpha_i成对优化的alpha_j                 j = selectJrand(i,m)                 #步骤1：计算误差Ej                 fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b                 Ej = fXj - float(labelMat[j])                 #保存更新前的aplpha值，使用深拷贝                 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[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                 #步骤3：计算eta                 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                 #步骤4：更新alpha_j                 alphas[j] -= labelMat[j]*(Ei - Ej)/eta                 #步骤5：修剪alpha_j                 alphas[j] = clipAlpha(alphas[j],H,L)                 if (abs(alphas[j] - alphaJold) < 0.00001):                     print("alpha_j变化太小");                     continue                 #步骤6：更新alpha_i                 alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])                 #步骤7：更新b_1和b_2                 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                 #步骤8：根据b_1和b_2更新b                 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("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))         #更新迭代次数         if (alphaPairsChanged == 0): iter_num += 1         else: iter_num = 0         print("迭代次数: %d" % iter_num)     return b,alphas   def showClassifer(dataMat, w, b):     '''     分类结果可视化     :param dataMat: 数据矩阵     :param w: 直线法向量     :param b: 直线解决     :return:     '''     #绘制样本点     data_plus = []                                  #正样本     data_minus = []                                 #负样本     for i in range(len(dataMat)):         if labelMat[i] > 0:             data_plus.append(dataMat[i])         else:             data_minus.append(dataMat[i])     data_plus_np = np.array(data_plus)              #转换为numpy矩阵     data_minus_np = np.array(data_minus)            #转换为numpy矩阵     plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)   #正样本散点图     plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图     #绘制直线     x1 = max(dataMat)[0]     x2 = min(dataMat)[0]     a1, a2 = w     b = float(b)     a1 = float(a1[0])     a2 = float(a2[0])     y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2     plt.plot([x1, x2], [y1, y2])     #找出支持向量点     for i, alpha in enumerate(alphas):         if abs(alpha) > 0:             x, y = dataMat[i]             plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')     plt.show()     def get_w(dataMat, labelMat, alphas):     '''     函数说明:计算w     :param dataMat: 数据矩阵     :param labelMat: 数据标签     :param alphas: alphas值     :return:     '''     alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)     w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)     return w.tolist()     if __name__ == '__main__':     dataMat, labelMat = loadDataSet('testSet.txt')     b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)     w = get_w(dataMat, labelMat, alphas)     showClassifer(dataMat, w, b)

网站代码和数据集已经上传到网站中，如果是微信请点击阅读原文下载，如果是知乎请转到个人网站在附件中下载。

我这里对代码做加强版说明：

其中容错率指的的是误差（即预测的结果和真实的结果进行对比，如果差别太大就对alpha值进行优化）不管是正间隔和负间隔都会被调整计算。

运行后如下：

我们发现实际上，上面的代码运行还是很消耗时间的，因为我们选择alpha的值是随机的，而我们期望的是对符合条件的误差项最大的alpha进行调整，我们如果用一个数据结构缓存每次更新的误差，那么速度将会大大的提高，具体代码在机器学习实战中有这段代码，这里就贴出来啦，后面的博客将会提到线性不可分，及其核函数的原理。