经过前几篇文章的学习,SVM的优化目标,SMO算法的基本实现步骤,模型对应参数的选择,我们已经都有了一定的理解,结合《机器学习实战》,动手实践一个基本的SVM支持向量机,来完成一个简单的二分类任务。
建立模型之前,首先看一下我们的数据,然后再用支持向量机实现分类:
这里只截取部分,我们的数据是一个txt文件,包含100个数据点,前两行为x1,x2,负一行为标签值,取值{-1,+1}.
读取文件数据
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#数据集,标签值
LoadDataSet函数的作用是读取数据文件,return得到两个列表,一个包含二维数据点,一个包含点对应的标签值
随机选择j
import random
def selectJrand(i,m):
j=i
while (j==i):
j = int(random.uniform(0,m))
return j
selectJrand函数通过while函数判断 j==i,从而达到随机生成不等于 i 的 j ,因为这里是简化的SVM算法,如果在复杂情况下,j的选择会和这里j所的选择方法有所不同
调整alpha
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
clipAlpha函数根据(4)的alpha的界限,调整迭代得到的alpha
简化版SMO算法
import numpy as np
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
#参数分别为数据集,类别标签,常数c,容错率和最大循环次数
#外循环
dataMatrix = np.mat(dataMatIn);labelMat = np.mat(classLabels).transpose() #数据与标签矩阵转化
b = 0;m,n = np.shape(dataMatrix)#读取数据维度
alphas = np.mat(np.zeros((m,1)))#初始化alpha为0
iter = 0#初始化迭代数
while (iter < maxIter):
#内循环
alphaParisChanged = 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)):#判断alpha是否满足KKT条件
j = selectJrand(i,m)#随机选取j
fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b #计算预测类别
Ej = fXj - float(labelMat[j])#计算误差e
alphaIold = alphas[i].copy()#为alphaIold和alphaJold分配新的内存
alphaJold = alphas[j].copy()
if (labelMat[i] != labelMat[j]):#alpha的上限与下限
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
#计算smo算法中的η
if eta > 0 : print("eta>=0");continue
#对i进行修改,修改量与j相同,但方向相反
alphas[j] -= labelMat[j]*(Ei-Ej)/eta
#通过L H的值优化
alphas[j] = clipAlpha(alphas[j],H,L)
#判断alphaJ的优化程度,更新alphaI,并根据新的alpha得到新的b
if (abs(alphas[j]-alphaJold) < 0.0001):
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
alphaParisChanged += 1
print("iter : %d i:%d,pairs changed %d" %(iter,i,alphaParisChanged))
if (alphaParisChanged == 0):#判断是否对alpha进行了更新
iter += 1
else:
iter = 0
print("iterration number :%d" % iter)
return b,alphas#返回超平面需要的参数alpha,b
SmoSimple函数是smo算法的简单实现版本,与原始算法思想基本一致,需要计算的是误差 ei = f(xi)-yi,二次项η (eta),以及alpha的上界H下界L,需要判断的是开头的KKT条件是否满足,判断yi与yj的类别是否一致从而决定上下界的形式,判断alpha的更新程度以及值的大小从而决定b,需要更新的是alpha的值,b的值,还有算法迭代次数,最后返回我们超平面构建需要的参数alpha,b . Tips: KKT条件的判断与西瓜书,《统计学习方法》有所不同,在(3)中解释了代码中的KKT条件是如何判断的.
绘制样本数据点,超平面与支持向量
import matplotlib.pyplot as plt
def plot_point(dataMat,labelMat,alphas,b):
for i in range(np.shape(dataMat)[0]):#绘制shape(dataMat)[0]个数据点
plt.scatter(dataMat[i][0],dataMat[i][1])
alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)#转为数组求w
sum = 0#初始化w
for i in range(np.shape(dataMat)[0]):
sum += alphas[i]*labelMat[i]*dataMat[i].T#通过求导后得到的w公式,求和得到w
# print(sum)
for i, alpha in enumerate(alphas):#根据kkt条件,选取alpha不为0的点作为支持向量
if abs(alpha) > 0:
x, y = dataMat[i]
plt.scatter(x, y, s=100, c = '', alpha=0.5, linewidth=1.5, edgecolor='red')
x = np.arange(0,10,0.01)
y = (sum[0]*x+b)/(-1*sum[1])#通过原始公式推导出超平面
plt.scatter(x,y,s=5,marker = 'h')
plt.show()
plot_point函数是负责绘制超平面的,书上只给了示意图,所以就只能自己动手丰衣足食了,首先绘制了所有数据点,由于数据点比较标准(无明显越界),所以为了好看就没有根据类别设置颜色,其次根据对偶问题得到w的计算公式
通过SmoSimple返回的alpha计算w,再加上b,就得到了最优分割超平面的公式,由于是二维数据点,所以稍加整理就可以得到超平面的二维表达式,然后我们设置适当的x范围,绘制超平面就ok了;还有一点就是支持向量的标注,因为alpha≠0时,(wi.T).x+b=0,对应的点在分割边界上,所以对应的点即为支持向量Support Vector,最后用c=' '空心标注出支持向量,就能得到书中的图了.
主函数
if __name__ == '__main__':
dataMat, labelMat = loadDataSet('testSet.txt')
b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
plot_point(dataMat,labelMat,alphas,b)
可以自己加print语句,看一下w,b:
b = [[-3.86378686]]
w = [ 0.79293619 -0.23660327]
[Finished in 12.9s]
然后看一下得到的分隔超平面:
总结:
第一次看书上的代码会有些懵,不过看一遍SMO算法推导,代码部分就基本无压力了,这里运用的是简化版算法,用时较长,不过得到的结果比较满意,之后还会用完整的SMO算法和Sklearn实现SVM,运行速度会大大提升,同时结果也比较满意。