目录
六、支持向量回归(Support Vector Regression,SVR)
一、最大间隔与分类
1.超平面介绍
我们可以看到下图的数据点分为两个类别。我们能够很容易地在途中画出一条直线将两组数据点分开,在这种情况下的,这组数据就被称为线性可分数据,而该直线就被称为超平面。
那么我们在往下看,可以看到有这么多条的直线都可以将两组数据分开,那么哪一条直线的效果最好呢,毫无疑问应该选择正中间红色那条,为什么呢?因为这条线相较于其他的线容忍性好、鲁棒性高、泛化能力最强,并且最大化了决策边界的边缘。
下图是在各个维度空间中SVM形式,我们可以看出如果一个空间维度是N维,那么在该空间上超平面的维度就应该是N-1。
2. 超平面方程
2.1 方程
2.2 方程的推导
假设二维空间下有一个直线方程为y=kx+b
变换该直线方程的形式为kx-y+b=0
然后将x和y替换成x1和x2,那么式子就变成了kx1-x2+b=0
再向量化得到列向量w=[a -1]和列向量x=[x1 x2]
那么式子最后就变成了上面的超平面方程。
2.3 最大间隔
我们在高中学过点到直线的距离,公式为
那么点到超平面的距离也是一样
那么分类间隔width=2d
在该副图中width = 2/||w||
2.4 分类器
输入数据给分类器会输出一个类别标签,利用单位阶跃函数,当输入数据大于0输出+1,小于0则输出-1.
假设我们的训练数据为:
线性可分当且仅当得出:
我们对上述式子进行变换,把yi乘上左边的超平面方程可以得到
接下来我们引入一个大于0的值c,因为上述式子大于0,c也大于0,所以上述式子可以大于或者等于c,于是得到
因为c是大于0的两边同时除以c,并且因为左边这个式子是一个很大的正数所以没什么影响,所以最后式子就变为
3. 最大化间隔
3.1 概述
我们需要找出分类器中定义的w和b,为此我们要找到具有最小间隔的数据点,找到之后对该间隔最大化,就可以写作:
对上述的式子直接求解非常困难,所以我们要将它转换成另一件更加容易求解的形式,如果令上述式子中min括号中的乘法都为1的话,那么就可以变成求||w||最小值来求最终解,但是事实并非如此,只有那些离超平面近的才为1,离超平面越远这个值也就越大。
我们要求解d的最大化问题转化乘求||w||的最小化问题,也就是
3.2 例子
给出三个点,求w1、w2、b的值。
二、对偶问题
1. 等式约束
给定一个目标函数 f : Rn→R,希望找到x∈Rn ,在满足约束条件g(x)=0的前提 下,使得f(x)有最小值。该约束优化问题记为:
可建立拉格朗日函数:L(x,λ) = f(x) + λg(x)
λ为拉格朗日乘数,将原本约束优化问题转换成等价的无约束优化问题:L(x,λ)
然后分别对x和λ求偏导,可以得到:
联立方程组就可以得出相应的解。
2. 不等式约束的KKT条件
将约束等式g(x)=0推广为不等式g(x)≤0,那么这个约束优化问题就变为:
定义可行域K,x*为满足约束条件的最佳解,那么分开成两种情况讨论:
1.g(x*)<0,最佳解位于K的内部,称为内部解,此时的约束条件是无效的。
内部解:因为约束条件无效,意思就是g(x)没有影响,所以驻点x*满足▽f=0且λ=0.
2.g(x*)=0,最佳解位于K的边界上,称为边界解,此时的约束条件是有效的。
边界解:约束条件有效,g(x)=0,在对式子求偏导我们得到▽f(x)+λ▽g(x) = 0,我们要求的是f(x)的最小值,因为g(x)≤0所以▽g(x)是指向K的外部的,所以我们可以得知λ≥0,这称为对偶可行性。
不论是哪种情况,λg(x)=0是恒成立的,这称为互补松弛性。
由上面我们就可以得出最佳解的必要条件为:
这些条件就被称为KKT条件。
3. 最大间隔问题的拉格朗日乘法
3.1 引入拉格朗日乘子构建拉格朗日函数
拉格朗日乘子αi≥0
3.2 分别对w和b求偏导,并令偏导为0
得到两个式子
3.3 优化后的式子
将w和b代回原式子可以得到
4. 支持向量机基本型
从最开始的求解w的最小值变成了求解αi的最大值
三、SMO算法
1. Platt的SMO算法
Platt的SMO算法是将大优化问题分解为多个小优化问题来求解的,这些小优化问题往往更容易求解,并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。因此它能做到在结果完全相同的同时,缩短计算的时间。
基本思路是不断执行如下两个步骤直至收敛:
1.选取一对需要更新的变量αi和αj。
2.固定除了αi和αj以外的参数,求解对偶问题更新αi和αj。
实际上SMO算法的工作原理就是每一次循环选择两个α进行优化处理,只要找到两个满足条件的α,就增大其中一个并同时减小另外一个,为什么要这样做呢?因为我们不能破环约束条件,只改变一个α可能会导致这个约束条件失效。这边的两个条件分别是1.这一对α要在间隔边界之外2.这一对α还没有进行过区间化处理或者不在边界上。
流程图如下:
假设最优解为
利用SMO算法我们就能得到w*和b*的值:
最后就得到我们想要求超平面方程
2. 简化版的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
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
b = 0; m,n = np.shape(dataMatrix)
alphas = np.mat(np.zeros((m,1)))
iter_num = 0
while (iter_num < maxIter):
alphaPairsChanged = 0
for i in range(m):
fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
Ei = fXi - float(labelMat[i])
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m)
fXj = float(np.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])
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_num,i,alphaPairsChanged))
if (alphaPairsChanged == 0): iter_num += 1
else: iter_num = 0
print("iteration number: %d" % iter_num)
return b,alphas
dataArr, labelArr = loadDataSet('svm_data/data/testSet.txt')
b,alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
输出支持向量。
for i in range(100):
if alphas[i]>0.0:
print(dataArr[i],labelArr[i])
3. 完整版的SMO算法
import matplotlib.pyplot as plt
import numpy as np
import random
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler, kTup):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2)))
self.K = np.mat(np.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(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def kernelTrans(X, A, kTup):
m,n = np.shape(X)
K = np.mat(np.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 = np.exp(K/(-1*kTup[1]**2))
else: raise NameError('Houston We Have a Problem That Kernal is not recognized')
return K
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 selectJ(i, oS, Ei):
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = np.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):
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):
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
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)
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.K[i,j] - oS.K[i,i] - oS.K[j,j]
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)
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])
updateEk(oS, i)
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]) 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, kTup = ('lin',0)):
oS = optStruct(np.mat(dataMatIn), np.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, pair changed %d" % (iter,i,alphaPairsChanged))
iter += 1
else:
nonBoundIs = np.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, pair 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
dataArr, labelArr = loadDataSet('svm_data/data/testSet.txt')
b,alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
四、核函数
1.线性不可分
如果在当前空间中找不到一个超平面来划分两类样本,应该怎么办呢?
这时候就出现了线性不可分的这种情况,之前我们讨论的都是线性可分的情况。这时就需要核函数,将样本从原始空间映射到一个更高维的特征空间,使得样本在这个特征空间内线性可分。
2.核支持向量机
假设样本x映射后的向量为Φ(x),划分超平面为
原始问题为
将上面公式用内积表示:
使用非线性映射,将数据映射到特征空间,也就是把上面式子的替换成,在特征空间使用线性学习器,变形完以后函数如下:
基本想法:不显式地构造核映射,而是设计核函数,实现。
Mercer定理(充分非必要):只要对称函数值所对应的核矩阵半正定,则该函数可作为核函数。
本来线性不可分的用核函数处理完过后,仍然可用SMO算法求解。
3. 常见的核函数
五、软间隔与正则化
1.软间隔的引入
在实际的应用中,很难选择合适的核函数,使样本在特征空间中线性可分,并且线性可分的结果也很难断定是否是由过拟合造成的。
由此我们引入软间隔,允许SVM在某些样本上不满足约束,意思也就是允许向量机在一些样本上出错。
部分的样本允许:
称为松弛变量
2.损失函数
其中C>0为惩罚参数;是“0/1损失函数”,当C无穷大时,会迫使所有样本满足约束。
因为非凸、非连续,所以常用其他的函数,称为替代损失,替代损失函数通常是凸的连续函数并且是的上界。
hinge loss:
exponenetail loss:
logistic loss:
3. 采用hinge loss作为替代损失函数
那么原始问题就表示为
引入松弛变量C,对偶问题就变为
上述优化目标满足KKT条件,可以推出最终的模型仅仅与支持向量有关系,说明hinge loss损失函数依然保持了支持向量机解的稀疏性。
4.一般形式
结构风险,用来描述模型的某一些性质
而后面是经验风险,用来描述模型与训练数据的契合程度
六、支持向量回归(Support Vector Regression,SVR)
1. 特点
允许模型输出和实际输出间存在2的误差,也就是如果训练样本落在该间隔带中会被认为预测正确。
2.损失函数
SVR问题可形式化为:
其中C是一个长度,为不敏感损失函数。
3. 形式化
七、实战:基于sklearn的垃圾邮件分类
利用sklearn库中svm来实现,首先分别读取数据集和测试集的数据,把特征值和标签值取出来,设置9个误差惩罚系数候选值,循环九次,判断这几个候选值哪一个更好,最后在用最好的一个C去构建模型,并且用模型预测输出的结果与测试集中的标签进行对比,输出每一次的真实值和预测值,并统计错误率。
import numpy as np
import pandas as pd
import scipy.io as sio
import matplotlib.pyplot as plt
from sklearn.svm import SVC
train_data = sio.loadmat('svm_data\\data\\spamTrain.mat')
test_data = sio.loadmat('svm_data\\data\\spamTest.mat')
X,y=train_data['X'],train_data['y']
Xtest,ytest=test_data['Xtest'],test_data['ytest']
Cvalues = [3, 10, 30, 100, 0.01, 0.03, 0.1, 0.3, 1]
best_score = 0
best_param = 0
for c in Cvalues:
svc = SVC(C=c, kernel='linear')
svc.fit(X, y.flatten())
score = svc.score(Xtest, ytest.flatten())
if score > best_score:
best_score = score
best_param = c
svc = SVC(C= best_param, kernel='linear')
svc.fit(X, y.flatten())
svc_predict = svc.predict(Xtest)
errorCount = 0
for i in range(len(Xtest)):
print("分类返回结果为%d\t真实结果为%d" % (ytest[i],svc_predict[i]))
if(ytest[i] != svc_predict[i]):
errorCount += 1.0
print("总共错了%d个数据\n错误率为%f%%" % (errorCount, errorCount/len(Xtest) * 100))
八、总结
在机器学习中,支持向量机(SVM)是一种很常用的分类算法,它主要是通过寻找一个最优的超平面来将不同类别的数据分开。当有许多个超平面都能分隔数据集时,应采用最大间隔,选择最佳的超平面,在这基础上可以引入SMO算法来计算w和b进而从w最小值变成最大化α。
在选择数据集的时候,数据集分为两种,一种是线性可分的数据集,另一种是线性不可分的数据集,线性不可分的数据集需要用到核函数。
构建完模型后,可以改变惩罚参数C来调整模型以达到最优的结果。
通过实验可以帮助我们更好地理解SVM模型的原理和性能,以及如何选择合适的数据集和参数来训练模型。这些实验结果可以帮助我们更好地应用SVM算法解决实际的分类问题。