Svm（support Vector Mac）又称为支持向量机，是一种二分类的模型。当然如果进行修改之后也是可以用于多类别问题的分类。支持向量机可以分为线性核非线性两大类。其主要思想为找到空间中的一个更够将所有数据样本划开的超平面，并且使得本本集中所有数据到这个超平面的距离最短。
一、基于最大间隔分隔数据
1.1支持向量与超平面
在了解svm算法之前，我们首先需要了解一下线性分类器这个概念。比如给定一系列的数据样本，每个样本都有对应的一个标签。为了使得描述更加直观，我们采用二维平面进行解释，高维空间原理也是一样。举个简单子：如下图所示是一个二维平面，平面上有两类不同的数据，分别用圆圈和方块表示。我们可以很简单地找到一条直线使得两类数据正好能够完全分开。但是能将据点完全划开直线不止一条，那么在如此众多的直线中我们应该选择哪一条呢？从直观感觉上看图中的几条直线，
图 1 图2
1.2寻找最大间隔
1.2.1点到超平面的距离公式
既然这样的直线是存在的，那么我们怎样寻找出这样的直线呢？与二维空间类似，超平面的方程也可以写成一下形式：
有了超平面的表达式之后之后，我们就可以计算样本点到平面的距离了。假设
其中||W||为超平面的范数，常数b类似于直线方程中的截距。
上面的公式可以利用解析几何或高中平面几何知识进行推导，这里不做进一步解释。
1.2.2最大间隔的优化模型
现在我们已经知道了如何去求数据点到超平面的距离，在超平面确定的情况下，我们就能够找出所有支持向量，然后计算出间隔margin。每一个超平面都对应着一个margin，我们的目标就是找出所有margin中最大的那个值对应的超平面。因此用数学语言描述就是确定w、b使得margin最大。这是一个优化问题其目标函数可以写成：
其中
为了后面计算的方便，我们将目标函数等价替换为：
这是一个有约束条件的优化问题，通常我们可以用拉格朗日乘子法来求解。拉格朗日乘子法的介绍可以参考这篇博客。应用拉格朗日乘子法如下：
令 
求L关于求偏导数得：
将(1.7)代入到(1.6)中化简得：
原问题的对偶问题为：
该对偶问题的KKT条件为
到此，似乎问题就能够完美地解决了。但是这里有个假设：数据必须是百分之百可分的。但是实际中的数据几乎都不那么“干净”，或多或少都会存在一些噪点。为此下面我们将引入了松弛变量来解决这种问题。
1.2.3松弛变量
由上一节的分析我们知道实际中很多样本数据都不能够用一个超平面把数据完全分开。如果数据集中存在噪点的话，那么在求超平的时候就会出现很大问题。从图三中课看出其中一个蓝点偏差太大，如果把它作为支持向量的话所求出来的margin就会比不算入它时要小得多。更糟糕的情况是如果这个蓝点落在了红点之间那么就找不出超平面了。
图 3
因此引入一个松弛变量ξ来允许一些数据可以处于分隔面错误的一侧。这时新的约束条件变为:
式中ξi的含义为允许第i个数据点允许偏离的间隔。如果让ξ任意大的话，那么任意的超平面都是符合条件的了。所以在原有目标的基础之上，我们也尽可能的让ξ的总量也尽可能地小。所以新的目标函数变为：
其中的C是用于控制“最大化间隔”和“保证大部分的点的函数间隔都小于1”这两个目标的权重。将上述模型完整的写下来就是：
新的拉格朗日函数变为：
接下来将拉格朗日函数转化为其对偶函数，首先对
代入原式化简之后得到和原来一样的目标函数：
但是由于我们得到
经过添加松弛变量的方法，我们现在能够解决数据更加混乱的问题。通过修改参数C，我们可以得到不同的结果而C的大小到底取多少比较合适，需要根据实际问题进行调节。
1.2.4核函数
以上讨论的都是在线性可分情况进行讨论的，但是实际问题中给出的数据并不是都是线性可分的，比如有些数据可能是如图4样子。
图4
那么这种非线性可分的数据是否就不能用svm算法来求解呢？答案是否定的。事实上，对于低维平面内不可分的数据，放在一个高维空间中去就有可能变得可分。以二维平面的数据为例，我们可以通过找到一个映射将二维平面的点放到三维平面之中。理论上任意的数据样本都能够找到一个合适的映射使得这些在低维空间不能划分的样本到高维空间中之后能够线性可分。我们再来看一下之前的目标函数：
定义一个映射使得将所有映射到更高维空间之后等价于求解上述问题的对偶问题：
这样对于线性不可分的问题就解决了，现在只需要找出一个合适的映射即可。当特征变量非常多的时候在，高维空间中计算内积的运算量是非常庞大的。考虑到我们的目的并不是为找到这样一个映射而是为了计算其在高维空间的内积，因此如果我们能够找到计算高维空间下内积的公式，那么就能够避免这样庞大的计算量，我们的问题也就解决了。实际上这就是我们要找的核函数
通过以上例子，我们可以很明显地看到核函数是怎样运作的。上述问题的对偶问题可以写成如下形式：
那么怎样的函数才可以作为核函数呢？下面的一个定理可以帮助我们判断。
Mercer定理：任何半正定的函数都可以作为核函数。其中所谓半正定函数
值得注意的是，上述定理中所给出的条件是充分条件而非充要条件。因为有些非正定函数也可以作为核函数。
下面是一些常用的核函数：
表1 常用核函数表
核函数名称 | 核函数表达式 | 核函数名称 | 核函数表达式 |
线性核 | 指数核 | ||
多项式核 | 拉普拉斯核 | ||
高斯核 | Sigmoid核 |
现在我们已经了解了一些支持向量机的理论基础，我们通过对偶问题的的转化将最开始求
二、Smo算法原理
2.1 约束条件
根据以上问题的分析，我们已经将原始问题转化为了求的值，即求下面优化模型的解：
求解
Smo算法的原理为：每次任意抽取两个乘子
而原对偶问题的子问题的目标函数可以表达成：
其中：
这里之所以算两个
为了表述方便，定义一个特征到输出结果的输出函数：
该对偶问题中KKT条件为：
根据上述问题的KKT条件可以得出目标函数中的的
1、 
2、
3、
最优解需要满足KKT条件，因此需要满足以上的三个条件都满足。而不满足这三个条件的情况也有三种：
1、
2、
3、
也就是说如果存在不满足这些KKT条件的，我们就要更新它，这就是约束条件之一。其次，
其中
2.2参数优化
因为两个因子不好同时求解，所以可先求第二个乘子
接下来，综合
当
当
回顾第二个约束条件 ：
令其两边同时乘y1，可得：
其中：
因此
令 
经过转化之后可得：
那么如何来选择乘子
而b在满足以下条件时需要更新：
且更新后的和如下：
每次更新完两个乘子之后，都需要重新计算b以及对应的E。最后更新完所有的
三、Svm的python实现与应用
第二节已经对Smo算法进行了充分的说明并进行了推导，现在一切都准备好了。接下来需要做的就是实现这些算法了，这里我参考了《机器学习实战》这本书中的代码，并利用该程序对一个问题进行了求解。由于代码数量过大，因此这里就不再列出，而是放在附录中。有兴趣的朋友可以直接下载，也可以去官网下载源代码。笔者在读这些代码的过程中，也遇到了许多困难，大部分都根据自己的情况进行了注释。
3.1问题描述
这里我选取的一个数据集为声呐数据集，该问题为需要根据声呐从不同角度返回的声音强度来预测被测物体是岩石还是矿井。数据集中共有208个数据，60个输入变量和1个输出变量。每个数据的前60个元素为不同方向返回的声音强度，最后一个元素为本次用于声呐测试的物体类型。其中标签M表示矿井，标签R为岩石。
3.2数据预处理
所给数据中没有缺失值和异常值，因此不需要对数据集进行清洗。注意到所给数据集的标签为字母类型，而svm算法的标准标签为“-1”和“1”的形式，所以需要对标签进行转化，用“-1”、“1”分别代替M和R。由于该数据集中所给标签相同的数据都放在了一起，因此先对数据顺序进行了打乱。代码如下：
def loadDataSet(fileName): dataMat=[];labelMat=[];data=[] fr=open(fileName) for line in fr.readlines(): line=line.replace(',','\t')#去除逗号 line=line.replace('M','-1')#对标签进行替换 line=line.replace('R','1') lineArr=line.strip('\n').split('\t')#分割数据 data.append([float(lineArr[i]) for i in range(len(lineArr))]) random.shuffle(data) #随机打乱列表 for i in range(len(data)): dataMat.append(data[i][0:len(data)-1]) labelMat.append(data[i][-1]) return dataMat,labelMat |
3.3模型训练及测试
首先测试了一下将所有数据即作为训练集又作为测试集，然后用Smo模型进行训练找到所有的支持向量。最后根据训练好的模型进行求解，最终测试的准确率平均为54%。如果把数据集分成测试集和训练集，发现测试的准确率也在这附近。而根据网上数据统计该数据集测试的准确率最高为84%，一般情况下为百分之六十几。数据集本身是造成测试准确率低的一个原因，但是另外一个更加重要的原因大概是参数的选择问题。如何选取合适的参数是一个值得思考的问题，在接下来的学习过程中我也会注意一下参数选取这个问题。到这里，关于svm的算法就告一段落了。
#svm算法的实现
from numpy import*
import random
from time import*
def loadDataSet(fileName):#输出dataArr(m*n),labelArr(1*m)其中m为数据集的个数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不同的jj=iwhile(j==i):j=int(random.uniform(0,m))return jdef clipAlpha(aj,H,L):#调整大于H或小于L的alpha的值if aj>H:aj=Hif aj<L:aj=Lreturn ajdef smoSimple(dataMatIn,classLabels,C,toler,maxIter):dataMatrix=mat(dataMatIn);labelMat=mat(classLabels).transpose()#转置b=0;m,n=shape(dataMatrix)#m为输入数据的个数，n为输入向量的维数alpha=mat(zeros((m,1)))#初始化参数，确定m个alphaiter=0#用于计算迭代次数while (iter<maxIter):#当迭代次数小于最大迭代次数时（外循环）alphaPairsChanged=0#初始化alpha的改变量为0for i in range(m):#内循环fXi=float(multiply(alpha,labelMat).T*\(dataMatrix*dataMatrix[i,:].T))+b#计算f(xi)Ei=fXi-float(labelMat[i])#计算f(xi)与标签之间的误差if ((labelMat[i]*Ei<-toler)and(alpha[i]<C))or\((labelMat[i]*Ei>toler)and(alpha[i]>0)):#如果可以进行优化j=selectJrand(i,m)#随机选择一个j与i配对fXj=float(multiply(alpha,labelMat).T*\(dataMatrix*dataMatrix[j,:].T))+b#计算f(xj)Ej=fXj-float(labelMat[j])#计算j的误差alphaIold=alpha[i].copy()#保存原来的alpha(i)alphaJold=alpha[j].copy()if(labelMat[i]!=labelMat[j]):#保证alpha在0到c之间L=max(0,alpha[j]-alpha[i])H=min(C,C+alpha[j]-alpha[i])else:L=max(0,alpha[j]+alpha[i]-C)H=min(C,alpha[j]+alpha[i])if L==H:print('L=H');continueeta=2*dataMatrix[i,:]*dataMatrix[j,:].T-\dataMatrix[i,:]*dataMatrix[i,:].T-\dataMatrix[j,:]*dataMatrix[j,:].Tif eta>=0:print('eta=0');continuealpha[j]-=labelMat[j]*(Ei-Ej)/etaalpha[j]=clipAlpha(alpha[j],H,L)#调整大于H或小于L的alphaif (abs(alpha[j]-alphaJold)<0.0001):print('j not move enough');continuealpha[i]+=labelMat[j]*labelMat[i]*(alphaJold-alpha[j])b1=b-Ei-labelMat[i]*(alpha[i]-alphaIold)*\dataMatrix[i,:]*dataMatrix[i,:].T-\labelMat[j]*(alpha[j]-alphaJold)*\dataMatrix[i,:]*dataMatrix[j,:].T#设置bb2=b-Ej-labelMat[i]*(alpha[i]-alphaIold)*\dataMatrix[i,:]*dataMatrix[i,:].T-\labelMat[j]*(alpha[j]-alphaJold)*\dataMatrix[j,:]*dataMatrix[j,:].Tif (0<alpha[i])and(C>alpha[j]):b=b1elif(0<alpha[j])and(C>alpha[j]):b=b2else:b=(b1+b2)/2alphaPairsChanged+=1print('iter:%d i:%d,pairs changed%d'%(iter,i,alphaPairsChanged))if (alphaPairsChanged==0):iter+=1else:iter=0print('iteraction number:%d'%iter)return b,alpha
#定义径向基函数
def kernelTrans(X, A, kTup):#定义核转换函数（径向基函数）m,n = shape(X)K = mat(zeros((m,1)))if kTup[0]=='lin': K = X * A.T #线性核K为m*1的矩阵elif kTup[0]=='rbf':for j in range(m):deltaRow = X[j,:] - AK[j] = deltaRow*deltaRow.TK = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlabelse: raise NameError('Houston We Have a Problem -- \That Kernel is not recognized')return Kclass optStruct:def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parametersself.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))) #first column is valid flagself.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):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):maxK = -1; maxDeltaE = 0; Ej = 0oS.eCache[i] = [1,Ei]validEcacheList = nonzero(oS.eCache[:,0].A)[0]if (len(validEcacheList)) > 1:for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta Eif k == i: continue #don't calc for i, waste of timeEk = calcEk(oS, k)deltaE = abs(Ei - Ek)if (deltaE > maxDeltaE):maxK = k; maxDeltaE = deltaE; Ej = Ekreturn maxK, Ejelse: #in this case (first time around) we don't have any valid eCache valuesj = selectJrand(i, oS.m)Ej = calcEk(oS, j)return j, Ejdef updateEk(oS, k):#after any alpha has changed update the new value in the cacheEk = calcEk(oS, k)oS.eCache[k] = [1,Ek]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) #added this for the 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
#smoP函数用于计算超平的alpha,b
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #完整的Platter SMOoS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)iter = 0#计算循环的次数entireSet = True; alphaPairsChanged = 0while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):alphaPairsChanged = 0if entireSet: #go over allfor i in range(oS.m):alphaPairsChanged += innerL(i,oS)print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))iter += 1else:#go over non-bound (railed) alphasnonBoundIs = 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 = Trueprint("iteration number: %d" % iter)return oS.b,oS.alphas
#calcWs用于计算权重值w
def calcWs(alphas,dataArr,classLabels):#计算权重WX = mat(dataArr); labelMat = mat(classLabels).transpose()m,n = shape(X)w = zeros((n,1))for i in range(m):w += multiply(alphas[i]*labelMat[i],X[i,:].T)return w#值得注意的是测试准确与k1和C的取值有关。
def testRbf(k1=1.3):#给定输入参数K1#测试训练集上的准确率dataArr,labelArr = loadDataSet('testSetRBF.txt')#导入数据作为训练集b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 importantdatMat=mat(dataArr); labelMat = mat(labelArr).transpose()svInd=nonzero(alphas.A>0)[0]#找出alphas中大于0的元素的位置#此处需要说明一下alphas.A的含义sVs=datMat[svInd] #获取支持向量的矩阵，因为只要alpha中不等于0的元素都是支持向量labelSV = labelMat[svInd]#支持向量的标签print("there are %d Support Vectors" % shape(sVs)[0])#输出有多少个支持向量m,n = shape(datMat)#数据组的矩阵形状表示为有m个数据，数据维数为nerrorCount = 0#计算错误的个数for i in range(m):#开始分类，是函数的核心kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))#计算原数据集中各元素的核值predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b#计算预测结果y的值if sign(predict)!=sign(labelArr[i]): errorCount += 1#利用符号判断类别### sign（a）为符号函数：若a>0则输出1，若a<0则输出-1.###print("the training error rate is: %f" % (float(errorCount)/m))#2、测试测试集上的准确率dataArr,labelArr = loadDataSet('testSetRBF2.txt')errorCount = 0datMat=mat(dataArr)#labelMat = mat(labelArr).transpose()此处可以不用m,n = shape(datMat)for i in range(m):kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + bif sign(predict)!=sign(labelArr[i]): errorCount += 1print("the test error rate is: %f" % (float(errorCount)/m))
def main():t1=time()dataArr,labelArr=loadDataSet('testSet.txt')b,alphas=smoP(dataArr,labelArr,0.6,0.01,40)ws=calcWs(alphas,dataArr,labelArr)testRbf()t2=time()print("程序所用时间为%ss"%(t2-t1))if __name__=='__main__':main()后记
这是第一次写博客，其中难免会出错，因此希望大家能够批评指正。首先非常感谢网上的一些朋友，在理解svm这算法他们给了我很多启发，在公式推导中给了我很多参考的地方。本文主要参考的资料是《机器学习实战》和《惊呼！理解svm的三种境界》这篇博客。对于svm虽然学的时间不长，但是对它一直有种特别的感觉。第一次听说svm是在做一个验证码识别问题的时候，但那时候我使用的是KNN算法，尽管效果还不错，但是我一直希望能够用svm算法来完成这个题目。本来这次是打算把这个问题一起解决的，但是由于时间关系，没有来得及做。只能等下次有空闲的时候再来做这个问题了。
