机器学习笔记（Chapter 06 - 支持向量机）

支持向量机（Support Vector Machineds，SVM）是一个二类问题的分类器，实现方法多样，这里采用了序列最小优化（SMO）实现方法，并通过核函数拓展到非线性可分的SVM。

SVM和最大边缘超平面

• SVM的优缺点
  • 优点：泛化错误率低，计算开销不大撒，结果易解释
  • 缺点：对参数调节和核函数的选择敏感，原始分类器不加修改情况下仅适用于处理二类问题
  • 适用数值类型：数值型和标称型
• 最大边缘超平面：在二维平面上分布的二类数值点，如果可以通过一条直线将两组不同类的数据分开，则这组数据线性可分。在假设数据线性可分的前提下，将数据集分开的直线被称为分隔超平面，如果数据分布在三位平面，那么分隔超平面就是二维的。如果数据集分布在N维空间，则分隔超平面是N-1维。如果数据点离分隔超平面越远，则最后的预测结果就越好。因为决策边界边缘较小的分类器对模型的过分拟合更加敏感，从而在未知的样本上的泛化能力很差。
• 支持向量：离分隔超平面最近的那些点，支持向量机决策只依赖支持向量。
• 寻找最大间隔：用向量的形式W·X+b书写分隔超平面不需要考虑空间维度，其中向量W和常量b描述了所给数据的分隔超平面。因此SVM需要寻找使分隔超平面成为最大边缘超平面的W和b。

分隔超平面目标函数的优化

• SVM工作原理：与Logistic回归类似，使用一个类似海维赛德阶跃函数的函数对所给数据的W·X+b的结果判定分类，如果结果大于0则输出+1，否则输出-1。使用+1和-1而不使用1和0的作用在于，可以通过一个统一公式来表示间隔或者数据点到分隔超平面的距离。

• 函数间隔和几何间隔：点到分隔超平面的函数间隔为y*(wx+b)，其中y是函数的类别标签（+1或-1）；点到超平面的几何间隔为y*(wx+b)/||w||。SVM使用几何间隔定义数据点和超平面的距离，因为如果使用函数间隔，则随着w的放大，（wx+b）的值也随之不断增大，此时最优化（最大化）距离无法确定w。《机器学习实战》对SVM的原理介绍很粗略，并且直接给出了最终的可以解决线性不可分情况的公式。《机器学习实战》中SMO之前的部分在July的支持向量机通俗导论的第一层有比较清楚的介绍。

  下面部分《机器学习实战》没有讲，在《数据挖掘导论》的5.5节。

• 边缘公式的优化：要最大化最小间隔几何距离，考虑离决策边界最近的数据，如果数据在决策边界上方，则wx+b的结果是正值，在下方为负值，我们可以固定一个因子，调整另一个因子来优化最大值。因此我们设置一个约束条件y*(wx+b)>=1，这意味着所有的数据都在wx+b>=1和wx+b<=-1的范围内，距离超平面越远的店，其wx+b的绝对值就越大，只有支持向量才满足y(wx+b)=1的。我们选取两个数据点，一个在wx+b=1直线上，一个在wx+b=-1直线上，相减得到w(x1-x2)=2，注意w、x1和x2都是向量，所以d=x1-x2就代表着两点之间平行于超平面法线方向的距离。因此d=2/||w||。要让d最大，等价于让f(w)=||w||^2/2最小。因此，调整后的目标函数是f(w)，并且受到y(wx+b)>=1的约束。目标函数是二次的，w和b是线性的，因此该问题是凸优化问题（凸函数一阶可微，二阶导衡非负），此时可以引入拉格朗日算子，并且根据KKT条件将不等式约束改为等式约束y(wx+b)-1=0 ，变为最小化Lp = ||w||^2/2 - ∑(λ(y(wx+b)-1)，观察这个式子，我们限定λ>=0。其一阶导数为0，得到w=∑λyx，∑λy=0。将这两个条件代入拉格朗日算子的公式中，就得到书中的最后的目标函数。

• 不可分情况的处理：如果有少数数据噪声，需要引入正值的松弛变量ε，修改约束条件为y(wx+b)-(1-ε)>=0，假设直线wx+b=-1+ε经过数据点P，并且平行于决策边界，那么P到wx+b=-1的距离是ε/||w||。因此，ε提供了决策边界在训练样本P上的误差估计。同样，因为我们在决策边界上允许了一定的错误，可能导致误分许多的实例，所以对松弛变量很大的边界进行惩罚，修改后的目标函数为f(w) = ||w||^2 /2 +C(∑ε)^k，其中C和k是用户指定的参数，用于对误分的数据进行惩罚。假定k=1。这样修改后问题的拉格朗日函数多了一项-∑με，利用KKT条件约束，一阶导数为0，得到额外条件μ+λ =C，因此0<=λ<=C，配合∑λy = 0，这就是书中最终给出的约束公式。

SMO求解最优化问题

推荐JerryLead博客中的支持向量机（五）SMO算法。

• SMO算法的目标是求出一系列α和b，这里的α就是上面约束条件中的λ（拉格朗日乘子），因为参考的博客和书中都用α，下面也都用α。只要求出了α，根据w=∑αyx，就能够求出w。工作原理是每次循环选择两个alpha进行优化处理，一旦找到一对可以优化的α，就增大其中一个，同时减少另外一个。这两个α的选择方法决定了SMO的效率和正确率。

• SMO算法里的辅助函数

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):
	j = i
	while (j == i):
		j = int(random.uniform(0,m))
	return j
	
def clipAlpha(aj, H, L):
	if aj > H:
		aj = H
	if L > aj:
		aj = L
	return aj

• 《机器学习实战》书中先给了简化版的SMO算法，每次先选定一个α，然后随机选取另一个α。如果所有向量都没有被优化，就增加迭代次数，直到达到要求的迭代次数。书中给出平均速度14.5s。

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].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 = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print "eta>=0"; continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
                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
                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 "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas

• 启发式选择方法：每次选择α时，优先选择样本前面系数0<α<C的α作优化，因为在界上（α为0或C）的样例对应的α一般不会更改。这种启发式搜索方法是选择第一个α用的，只要选择出来的两个α中有一个违背了KKT条件，那么目标函数在一步迭代后值会减小。违背KKT条件不代表0<α<C，在界上也有可能会违背。因此在给定初始值α1=0后，先对所有样例进行循环，循环中碰到违背KKT条件的（不管界上还是界内）都进行迭代更新。等这轮过后，如果没有收敛，第二轮就只针对的样例进行迭代更新。在第一个α选择后，第二个α也使用启发式方法选择，第二个α的迭代步长大致正比于|E1-E2|，选择第二个α能够最大化|E1-E2|。即当E1为正时选择负的绝对值最大的E2，反之，选择正值最大的E2。最后的收敛条件是在界内（0<α<C）的样例都能够遵循KKT条件，且其对应的α只在极小的范围内变动。
• 完整的Platt SMO算法，书上数据平均时间0.78秒，下面是用到的辅助函数和结构。

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler):  # Initialize the structure with the parameters 
        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))) #first column is valid flag
        
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
        
def selectJK(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    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 E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEkK(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

• 完整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 selectJrand
        alphaIold = 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 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
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        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
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0

• 下面是外循环代码

def smoP(dataMatIn, classLabels, C, toler, maxIter):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)
                print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        else:#go over non-bound (railed) alphas
            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 #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print "iteration number: %d" % iter
    return oS.b,oS.alphas

• 下面是求W和分类函数

def calcWs(alphas, dataArr, classLabels): 
	X = mat(dataArr); labelMat = mat(classLabels).transpose()
	m, n = shape(X)
	w = zeros((n,1))	# n*1
	for i in range(m):
		w += multiply(alphas[i] * labelMat[i], X[i,:].T)	# n*1*1 (1*n)^T
	return w

def classified(dat, ws, b):
	if dat * mat(ws) + b > 0:
		return 1
	else:
		return -1

核函数

来自《数据挖掘导论》，并参考知乎上关于机器学习中核函数的讨论

• 径向基函数（RBF）：是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。

• 核函数和SVM是两个正交的概念，通过核函数可以将当前维度无法线性划分的数据转移到高维（无穷维度）。SVM核的变换后空间也称为再生核希尔伯特空间（RKHS），使用核函数计算点积开销更小，并且计算在原空间进行，无须担心维灾难问题。

• Mercer定理：对非线性SVM使用的核函数的主要要求是，必须存在一个相应的变换，使得计算一对向量的核函数等价于在变换后的空间中计算这对向量的点积。核函数K可以表示为K(u, v) = Φ(u)Φ(v)，当且仅当对于任意满足∫g(x)^2dx为有限值得函数g(x)，则∫K(x,y)g(x)g(y)dxdy >= 0。满足这个定理的核函数称为正定核函数。例如K(x,y) = (x·y+1)^p，K(x,y) = e^(-(|x-y|^2)/2σ^2))，K(x,y) = tanh(ky·y - δ)。
• 核函数转换

def kernelTrans(X, A, kTup):
	m, n = shape(X)
	K = mat(zeros((m,1)))
	if kTup[0] == 'lin': K = X*A.T
	elif kTup[0] == 'rbf':
		for j in range(m):
			deltaRow = X[j,:] - A
			K[j] = deltaRow * deltaRow.T
		K = exp(K / (-1*kTup[1]**2))
	else: raise NameError("Houston We Have a Problem -- \		That kernel is not recognized")
		That kernel is not recognized"
	return K

• 下面是测试函数，需要对函数innerL和calcEk和类optStruct做一定修改。

def testRbf(k1=1.3):
	dataArr, labelArr = loadDataSet('testSetRBF.txt')
	b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf',k1))
	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,:], ('rbf',k1))
		predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
		if sign(predict) != sign(labelArr[i]) : errorCount += 1
	print "the training error rate is %f" % (float(errorCount)/m) 
	dataArr, labelArr = loadDataSet('testSetRBF2.txt')
	datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
	m, n = shape(datMat)
	errorCount = 0
	for i in range(m):
		kernelEval = kernelTrans(sVs, datMat[i,:], ('rbf',k1))
		predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
		if sign(predict) != sign(labelArr[i]) : errorCount += 1
	print "the test error rate is %f" % (float(errorCount)/m)

• 原代码需要修改的地方

def innerL():
        ···
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
        ···
        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]
    	···
    	

def calcEk(oS, k):
    	fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:,k] + oS.b)
    	Ek = fXk - float(oS.labelMat[k])
    	return Ek
    	

class optStruct:
	def __init__(self, dataMatIn, classLabels, C, toler, kTup):
		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)))
		self.K = mat(zeros((self.m, self.m)))
		for i in range(self.m):
			self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

kNN手写问题回顾

• SVM是二类分类器，将非9的数字判为-1，否则判为1。
• Code - testDigits - svmMLiA.py

def img2vector(filename):
    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 returnVect

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)           #load the training set
    m = len(trainingFileList)
    trainingMat = zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     #take off .txt
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels

def testDigits(kTup=('rbf', 10)):
    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 += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadImages('testDigits')
    errorCount = 0
    datMat=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]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)

• 修改径向基核函数的参数σ，观察错误率。σ下降，则训练错误率降低，测试错误率上升。最小的训练错误率并不对应于最小的向量支持数目。另外线性和函数的效果并不很糟糕，可以牺牲线性核函数的错误率来换取分类速度的提高。

多类分类问题

• 第一种方法将多类问题分解为K个二类问题，对于类别yi，属于类别yi的为一类，不属于yi的为另一类。此方法称为一对其他（1-r）方法。
• 第二种方法为一对一（1-1）方法，构建K(K-1)/2个分类器，没一个分类器用来区分一对类(yi,yj)，此时忽略其他类的样本。
• 以上两种方法都使用二类分类器的组合预测，并投票表决，票数多的分类为最终分类。这种方法可能导致不同类的平局。
• 纠错输出编码：1-r和1-1方法都对二元分类的错误太敏感。可以参考海明编码，为每个类别分配一个码字，码字的每个二进制位训练一个二元分类器。

支持向量机总结

支持向量机是一种二类分类器，通过求解一个二次优化问题来最大化分类间隔。通过SMO算法每次优化两个α可以提升SVM的训练速度。核函数可以从一个低维空间的非线性数据映射到一个高维空间的线性数据i，此部分可以参考知乎。

SVM问题可以表示为凸优化问题，利用已知的有效算法发现目标函数的全局最小值。通过最大化决策边界的边缘来控制模型的能力，用户必须提供其他参数，如核函数类型、松弛变量带来的惩罚C。

---

参考文献： 《机器学习实战 - 美Peter Harrington》、《数据挖掘导论 - 美Pang-Ning Tan等》

参考的文章等：July的文章，JerryLead的文章，王赟 Maigo等在知乎上的答案

原创作品，允许转载，转载时无需告知，但请务必以超链接形式标明文章原始出处(https://forec.github.io/2016/02/11/machinelearning6/) 、作者信息（Forec）和本声明。