支持向量机 SVM

一、什么是支持向量机

支持向量机（support vector machines, SVM）是一种二分类模型，它的基本模型是定义在特征空间上的间隔最大的线性分类器，间隔最大使它有别于感知机；SVM还包括核技巧，这使它成为实质上的非线性分类器。SVM的的学习策略就是间隔最大化，可形式化为一个求解凸二次规划的问题，也等价于正则化的合页损失函数的最小化问题。SVM的的学习算法就是求解凸二次规划的最优化算法。

二、支持向量机算法原理

1.SVM基本思想

SVM学习的基本想法是求解能够正确划分训练数据集并且几何间隔最大的分离超*面。如下图所示，w*x+b=0 即为分离超*面，对于线性可分的数据集来说，这样的超*面有无穷多个（即感知机），但是几何间隔最大的分离超*面却是唯一的。

2.最大间隔与分类

SVM最大分类间隔的灵感来自于一个非常符合直觉的观察，如果存在两类数据，数据的特征是二维的，那么我们就可以把数据画在一个二维*面上，此时我想找到一个决策面（决策边界）去将这两类数据分开。如下图所示：

理论上这个决策边界有无数种选择，就像图中画出的四条黑色的线，都能实现分类，但是哪一种是最好的分类方式呢？SVM算法认为在上图中靠*决策*边界的点（正负样本）与决策边界的距离最大时，是最好的分类选择：

上图中红色的线就是要优化的目标，它表征了数据到决策边界的距离，这个距离就是所谓的最大分类间隔。同时在上面的几个数据，如果靠*两侧的数据少了几个，也不会影响决策边界的确定，而被红色框框出来三个数据才决定了最终的决策边界，所以这三个数据被称之为支持向量。

3.拉格朗日乘子法与对偶问题

很显然上述问题是一个凸二次规划问题，可以先对其使用拉格朗日乘子法得到其对偶问题，对每条约束添加拉格朗日乘子，由于这里的不是等式约束，所以每条约束\large \alpha_{i} 不同于等式的取值范围，后续会证明，如下式为该问题的拉格朗日函数：

分别对W和b求偏导数，得到：

4.SMO概念

SMO算法是一种启发式算法，其基本思路是：如果所有变量的解都满足此优化问题的KKT条件，那么这个最优化问题的解就得到了。(KKT条件是该最优化问题的充分必要条件)。否则，选择两个变量，固定其他变量针对这两个变量构建一个二次规划问题。

三、SVM代码实现

1.数据集的准备：

2.代码：

from numpy import *

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 #返回数据特征和数据类别

def selectJrand(i,m): #在0-m中随机选择一个不是i的整数
j=i
while (j==i):
j=int(random.uniform(0,m))
return j

def clipAlpha(aj,H,L): #保证a在L和H范围内（L <= a <= H）
if aj>H:
aj=H
if L>aj:
aj=L
return aj

def kernelTrans(X, A, kTup): #核函数，输入参数,X:支持向量的特征树；A：某一行特征数据；kTup：('lin',k1)核函数的类型和参数
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]'lin': #线性函数
K = X * A.T
elif kTup[0]'rbf': # 径向基函数(radial bias function)
for j in range(m):
deltaRow = X[j,:] - A
K[j] = deltaRowdeltaRow.T
K = exp(K/(-1kTup[1]**2)) #返回生成的结果
else:
raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
return K

定义类，方便存储数据

class optStruct:
def init(self,dataMatIn, classLabels, C, toler, kTup): # 存储各类参数
self.X = dataMatIn #数据特征
self.labelMat = classLabels #数据类别
self.C = C #软间隔参数C，参数越大，非线性拟合能力越强
self.tol = toler #停止阀值
self.m = shape(dataMatIn)[0] #数据行数
self.alphas = mat(zeros((self.m,1)))
self.b = 0 #初始设为0
self.eCache = mat(zeros((self.m,2))) #缓存
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): #计算Ek（参考《统计学习方法》p127公式7.105）
fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek

随机选取aj，并返回其E值

def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = nonzero(oS.eCache[:,0].A)[0] #返回矩阵中的非零位置的行数
if (len(validEcacheList)) > 1:
for k in validEcacheList:
if k == i:
continue
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE): #返回步长最大的aj
maxK = k
maxDeltaE = deltaE
Ej = Ek
return maxK, Ej
else:
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej

def updateEk(oS, k): #更新os数据
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]

首先检验ai是否满足KKT条件，如果不满足，随机选择aj进行优化，更新ai,aj,b值

def innerL(i, oS): #输入参数i和所有参数数据
Ei = calcEk(oS, i) #计算E值
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)): #检验这行数据是否符合KKT条件 参考《统计学习方法》p128公式7.111-113
j,Ej = selectJ(i, oS, Ei) #随机选取aj，并返回其E值
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
if (oS.labelMat[i] != oS.labelMat[j]): #以下代码的公式参考《统计学习方法》p126
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 LH:
print("LH")
return 0
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #参考《统计学习方法》p127公式7.107
if eta >= 0:
print("eta>=0")
return 0
oS.alphas[j] -= oS.labelMat[j](Ei - Ej)/eta #参考《统计学习方法》p127公式7.106
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) #参考《统计学习方法》p127公式7.108
updateEk(oS, j)
if (abs(oS.alphas[j] - alphaJold) < oS.tol): #alpha变化大小阀值（自己设定）
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]oS.labelMat[i](alphaJold - oS.alphas[j])#参考《统计学习方法》p127公式7.109
updateEk(oS, i) #更新数据
#以下求解b的过程，参考《统计学习方法》p129公式7.114-7.116
b1 = 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]<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:
return 0

SMO函数，用于快速求解出alpha

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #输入参数：数据特征，数据类别，参数C，阀值toler，最大迭代次数，核函数（默认线性核）
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
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)
print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) #显示第多少次迭代，那行特征数据使alpha发生了改变，这次改变了多少次alpha
iter += 1
else:
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 += 1
if entireSet:
entireSet = False
elif (alphaPairsChanged == 0):
entireSet = True
print("iteration number: %d" % iter)
return oS.b,oS.alphas

def testRbf(data_train,data_test):
dataArr,labelArr = loadDataSet(data_train) #读取训练数据
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 1.3)) #通过SMO算法得到b和alpha
datMat=mat(dataArr)
labelMat = mat(labelArr).transpose()
svInd=nonzero(alphas)[0] #选取不为0数据的行数（也就是支持向量）
sVs=datMat[svInd] #支持向量的特征数据
labelSV = labelMat[svInd] #支持向量的类别（1或-1）
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,:],('rbf', 1.3)) #将支持向量转化为核函数
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b #这一行的预测结果（代码来源于《统计学习方法》p133里面最后用于预测的公式）注意最后确定的分离*面只有那些支持向量决定。
if sign(predict)!=sign(labelArr[i]): #sign函数 -1 if x < 0, 0 if x==0, 1 if x > 0
errorCount += 1
print("the training error rate is: %f" % (float(errorCount)/m)) #打印出错误率
dataArr_test,labelArr_test = loadDataSet(data_test) #读取测试数据
errorCount_test = 0
datMat_test=mat(dataArr_test)
labelMat = mat(labelArr_test).transpose()
m,n = shape(datMat_test)
for i in range(m): #在测试数据上检验错误率
kernelEval = kernelTrans(sVs,datMat_test[i,:],('rbf', 1.3))
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
if sign(predict)!=sign(labelArr_test[i]):
errorCount_test += 1
print("the test error rate is: %f" % (float(errorCount_test)/m))

主程序

def main():
filename_traindata='C:\Users\Administrator\Desktop\data\traindata.txt'
filename_testdata='C:\Users\Administrator\Desktop\data\testdata.txt'
testRbf(filename_traindata,filename_testdata)

if name=='main':
main()

3.实验结果：

四、SVM的优缺点

优点：

1.支持向量机算法可以解决小样本情况下的机器学习问题，简化了通常的分类和回归等问题。
2.由于采用核函数方法克服了维数灾难和非线性可分的问题，所以向高维空间映射时没有增加计算的复杂性。换句话说，由于支持向量计算法的最终决策函数只由少数的支持向量所确定，所以计算的复杂性取决于支持向量的数目，而不是样本空间的维数。
3.支持向量机算法利用松弛变量可以允许一些点到分类*面的距离不满足原先要求，从而避免这些点对模型学习的影响。

缺点：

1.支持向量机算法对大规模训练样本难以实施。这是因为支持向量机算法借助二次规划求解支持向量，这其中会涉及m阶矩阵的计算，所以矩阵阶数很大时将耗费大量的机器内存和运算时间。
2.经典的支持向量机算法只给出了二分类的算法，而在数据挖掘的实际应用中，一般要解决多分类问题，但支持向量机对于多分类问题解决效果并不理想。
3.SVM算法效果与核函数的选择关系很大，往往需要尝试多种核函数，即使选择了效果比较好的高斯核函数，也要调参选择恰当的参数。另一方面就是现在常用的SVM理论都是使用固定惩罚系数C，但正负样本的两种错误造成的损失是不一样的。