机器学习实战之支持向量机（Support Vector Machine）

向量机（Support Vector Machine）"/>

机器学习实战之支持向量机（Support Vector Machine）

机器学习实战之支持向量机（Support Vector Machine）

理论知识

支持向量机分为：

• 线性可分支持向量机
• 线性支持向量机
• 非线性支持向量机

支持向量机的本质为求解一个凸二次规划问题，如今使用最广泛的求法为SMO（Sequential Minimal optimization）算法

在这里，不过多介绍理论知识，主要介绍实战内容

SVM应用的一般框架

1. 收集数据
2. 准备数据：数值型数据
3. 分析数据
4. 训练算法：SVM的大部分时间源自训练，该过程主要实现两个参数的调优
5. 测试算法
6. 使用算法：适用于几乎所有分类问题。

值得注意的是，SVM本身为一个二分类算法，如果要进行多分类，需要对SVM算法做一定修改。

SMO高效优化算法

SVM本质上是在求解一个凸二次规划问题，可以用常规的求解二次规划问题的算法来求解，但Platt提出的SMO是一个更高效的求解算法。

Platt的SMO算法

SMO算法是将大优化问题分解为多个小优化问题来求解。这些小优化问题往往易于求解，并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。（类似于分治算法）

其工作原理为：

每次循环中选择两个alpha进行优化处理。

一旦找到一对合适的alpha，就增大其中一个同时减少另一个。

应用简化版SMO算法处理小规模数据集

简化版SMO选择两个alpha的规则为：

首先在数据集上遍历每一个alpha，然后在剩下的alpha集合总随机选择另一个alpha，从而构成alpha对。

其伪代码为：

 创建一个alpha向量并将其初始化为0向量当迭代次数小于最大迭代次数时（外循环）对数据集中的每个数据向量（内循环）：如果该数据向量可以被优化：随机选择另一个数据向量同时优化这两个向量如果两个向量都不能被优化，退出内循环如果所有向量都没被优化，增加迭代次数，继续下一次循环。

一些辅助函数：

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, labelMatdef selectJrand(i, m):"""i为第一个alpha的下标m为所有alpha的数目即随机选择第二个alpha"""j = iwhile (j == i):j = int(random.uniform(0, m))return jdef clipAlpha(aj, H, L):"""用于调整大于H活小于L的alpha值"""if aj > H:aj = Hif L > aj:aj = Lreturn aj

简化版SMO算法代码如下：

from numpy import *def smoSimple(dataMatIn, classLabels, C, toler, maxIter):"""dataMatIn: 数据集classLabels: 类别标签Ctoler: 容错率maxIter: 退出前的最大迭代次数"""dataMatrix = mat(dataMatIn)labelMat = mat(classLabels).transpose() # 类别标签转换为列向量b = 0m, n = shape(dataMatrix)alphas = mat(zeros((m, 1)))  # alpha向量，存储所有alphaiter = 0  # 存储没有任何alpha改变的情况下遍历数据集的次数while (iter < maxIter):alphaPairsChanged = 0  # 记录alpha是否已经进行优化for i in range(m):fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b  # 预测的类别Ei = fXi - float(labelMat[i])  # 预测结果与真实结果的误差Ei# 如果误差较大，对alpha进行优化if ((labelMat[i] * Ei < -toler) and (alphas[i] < C) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0))):j = selectJrand(i, m)  # 随机选择第二个alpha值fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b  # Ej = fXj - float(labelMat[j])  # alphaIold = alphas[i].copy()alphaJold = alphas[i].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为alphas[j]的最佳修改量eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - dataMatrix[i, :] * dataMatrix[i, :].T - dataMatrix[j, :] * dataMatrix[j, :].Tif eta >= 0:print("eat >= 0")continuealphas[j] -= labelMat[j] * (Ei - Ej) / etaalphas[j] = clipAlpha(alphas[j], H, L) # 利用L和H对alphas[j]进行调整if (abs(alphas[j] - alphaJold) < 0.00001):print("j not moving enough")continuealphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])  # alphas[i]与alphas[j]的改变方向相反b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[i, :].T - \labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].Tb2 = b - Ej -labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - \labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].Tif (0 < alphas[i]) and (C > alphas[i]):b = b1elif (0 < alphas[j]) and (C > alphas[j]):b = b2else:b = (b1 + b2) / 2.0alphaPairsChanged += 1print("iter:%d i:%d, pairs change %d" % (iter, i, alphaPairsChanged))# 如果alpha对无更新，则迭代次数加一if (alphaPairsChanged == 0):iter += 1# 如果alpha对更新，则重新开始迭代else:iter = 0print("iteration number: %d" % iter)return b, alphas

利用完整版Platt SMO算法加速优化

简化版SMO算法在小样本数据集上效果还可以，但在大样本数据集上执行效率以及结果并不理想。

所以对简化版SMO算法加以改进，得到完整版的Platt SMO算法。其核心为利用一种启发式方法选取alpha

其基本原理为：

通过一个外循环来选择第一个alpha，其方式由两种：

1. 在所有数据集上进行单边扫描
2. 在边界alpha中进行单边扫描（非边界指的是不等于边界0或C的alpha值

确定一个alpha后，通过一个内循环来选择第二个alpha值。

​ 通过最大化步长的方式来获得第二个alpha值，选择使得步长或者 E i − E j E_i-E_j Ei​−Ej​最大的alpha值。

• 一些辅助函数：

class optStruct:def __init__(self, dataMatIn, classLabels, C, toler, kTup):"""使用一个类来封装算法所需的变量"""self.X = dataMatInself.labelMat = classLabelsself.C = Cself.tol = tolerself.m = shape(dataMatIn)[0]self.alphas = mat(zeros((self.m, 1)))self.b = 0self.eCache = mat(zeros((self.m, 2)))  # eCache为m * 2的矩阵，第一列为是否有效的标志位，第二列给出实际的E值self.K = mat(zeros((self.m, self.m)))for i in range(self.m):self.K[:, i] = kernelTrans(self.X, self.X[i,:], kTup)def calcEk(oS, k):"""计算E值并返回"""fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)Ek = fXk - float(oS.labelMat[k])return Ekdef selectJ(i,oS, Ei):"""选择第二个alpha的启发式方法"""maxK =  -1maxDeltaE = 0Ej = 0oS.eCache[i] = [1, Ei]validEcacheList = nonzero(oS.eCache[:, 0].A)[0] # 返回eCache第一列非零元素的索引值，从行的维度# .A 将矩阵转化为一个numpy数组，作为nonzero()的参数，因为矩阵无法作为nonzero()的参数if (len(validEcacheList)) > 1:# 选择使得E改变最大的 作为jfor k in validEcacheList:if k == i:continueEk = calcEk(oS, k)deltaE = abs(Ei - Ek)if (deltaE > maxDeltaE):maxK = kmaxDeltaE = deltaEEj = Ekreturn maxK, Ejelse:j = selectJrand(i, oS.m)Ej = calcEk(oS, j)return j, Ejdef updateEk(oS, k):"""计算误差值并存储到eCache中"""Ek = calcEk(oS, k)oS.eCache[k] = [1, Ek]

• Platt SMO的优化例程：

def innerL(i, oS):Ei = calcEk(oS, i)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)):j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrandalphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();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 0eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernelif eta >= 0: print ("eta>=0"); return 0oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/etaoS.alphas[j] = clipAlpha(oS.alphas[j],H,L)updateEk(oS, j) # alpha改变时更新eCacheif (abs(oS.alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); return 0oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as jupdateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie directionb1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2else: oS.b = (b1 + b2)/2.0return 1else: return 0

• 完整版的Platt SMO算法代码：

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMOoS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)iter = 0entireSet = True  # 标识是否遍历了整个集合alphaPairsChanged = 0  # 标识是否对alpha做了修改# 当迭代次数超过最大值，或者遍历整个集合都未对alpha做出修改，退出循环while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):alphaPairsChanged = 0if entireSet:  # 遍历数据集上任意可能的alphafor i in range(oS.m):        alphaPairsChanged += innerL(i,oS)print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))iter += 1else:# 遍历所有的非边界alpha值，即不再边界0或C上的值nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]for i in nonBoundIs:alphaPairsChanged += innerL(i,oS)print ("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))iter += 1if entireSet: entireSet = False #toggle entire set loopelif (alphaPairsChanged == 0): entireSet = True  return oS.b,oS.alphas

在复杂数据上应用核函数

上述的讨论都是建立在线性可分的条件下的，但现实中，样本数据往往是线性不可分的。这时候，我们就需要引入核函数来处理线性不可分情况。

核函数的作用：

将样本所在特征空间映射到另一个特征空间。

实现将低维特征空间的线性不可分情况转换为高维特征空间的线性可分情况，从而进行处理。

• 常用的核函数有：

  • 径向基核函数（又称高斯核函数）

  k ( x , y ) = e − ∣ ∣ x − y ∣ ∣ 2 2 2 σ 2 k(x,y)=e^{\frac{-||x-y||_2^2}{2\sigma^2}} k(x,y)=e2σ2−∣∣x−y∣∣22​​

• 核函数代码：

def kernelTrans(X, A, kTup):"""输入X A， 返回核函数的计算结果kTup为一个包含核函数信息的元组，第一个参数为描述所用核函数类型的一个字符串，其他2个参数则都是核函数可能需要的可选参数"""m, n = shape(X)K = mat(zeros((m, 1)))if kTup[0] == 'lin':K = X * A.Telif kTup[0] == 'rbf':for j in range(m):deltaRow = X[j, :] - AK[j] = deltaRow * deltaRow.TK = exp(K / (-1 * kTup[1] ** 2))  #else:raise NameError('Houston We have a problem -- That Kernel is not recogized')return K

在测试中使用核函数

使用径向基核函数处理线性不可分情况。

代码如下：

def testRbf(k1=1.3):dataArr, labelArr = loadDataSet('testSetRBF.txt')b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))dataMat = mat(dataArr)labelMat = mat(labelArr).transpose()sVInd = nonzero(alphas.A > 0)[0]  # 支持向量的索引sVs = dataMat[sVInd]  # 支持向量labelSV = labelMat[sVInd]  # 支持向量的类别标签print("there are %d Support Vectors" % shape(sVs)[0])m, n = shape(dataMat)errorCount = 0# 训练集for i in range(m):kernelEval = kernelTrans(sVs, dataMat[i, :], ('rbf', k1))predict = kernelEval.T * multiply(labelSV, alphas[sVInd]) + b  # 预测结果if sign(predict) != sign(labelArr[i]):errorCount += 1print("the training error rate is: %f" % (float(errorCount) / m))dataArr,labelArr = loadDataSet('testSetRBF2.txt')errorCount = 0dataMat=mat(dataArr);labelMat = mat(labelArr).transpose()m,n = shape(dataMat)for i in range(m):kernelEval = kernelTrans(sVs,dataMat[i,:],('rbf', k1))predict=kernelEval.T * multiply(labelSV,alphas[sVInd]) + bif sign(predict)!=sign(labelArr[i]): errorCount += 1 print("the test error rate is: %f" % (float(errorCount) / m))

示例：手写识别问题回顾

其过程为：

1. 收集数据：图片的文本文件
2. 准备数据：基于二值图像构造向量
3. 分析数据：
4. 训练算法：采用两种不同的核函数，并对径向基核函数采用不同的设置来运行SMO算法
5. 测试算法：编写一个函数来测试不同的核函数并计算错误率
6. 使用算法

其代码如下：

def img2vector(filename):# 将表示图片的(32,32)矩阵转换为(1, 1024)returnVect = zeros((1, 1024))fr = open(filename)for i in range(32):lineStr = fr.readline()for j in range(32):returnVect[0, 32 * i + j] = int(lineStr[j])return returnVectdef loadImages(dirName):"""贴标签，9为-1，其他为1"""from os import listdirhwLabels = []trainingFileList = listdir(dirName)m = len(trainingFileList)trainingMat = zeros((m, 1024))for i in range(m):fileNameStr = trainingFileList[i]fileStr = fileNameStr.split('.')[0]classNumStr = int(fileStr.split('-')[0])if classNumStr == 9:hwLabels.append(-1)else:hwLabels.append(1)trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))return trainingMat, hwLabelsdef testDigits(kTup=('rbf', 10)):"""测试使用径向基核函数的SVM算法"""dataArr, labelArr = loadImages('trainingDigits')b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)datMat = mat(dataArr)labelMat = mat(labelArr).transpose()svInd = nonzero(alphas.A > 0)[0]sVs = datMat[svInd]labelSV = labelMat[svInd]print("there are %d Support Vectors" % shape(sVs)[0])m, n = shape(datMat)errorCount = 0# 训练集for i in range(m):kernelEval = kernelTrans(sVs, datMat[i, :], kTup)predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b  # 预测结果if sign(predict) != sign(labelArr[i]):errorCount += 1print("the training error rate is: %f" % (float(errorCount) / m))dataArr,labelArr = loadImages('testDigits')errorCount = 0datMat=mat(dataArr);labelMat = mat(labelArr).transpose()m,n = shape(datMat)for i in range(m):kernelEval = kernelTrans(sVs,datMat[i,:],kTup)predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + bif sign(predict)!=sign(labelArr[i]): errorCount += 1 print("the test error rate is: %f" % (float(errorCount) / m))