SVM有很多实现,本篇关注其中最常用的一种,即序列最小优化(Sequential Minimal Optimization, 即SMO,算法的数学逻辑上一篇有介绍)算法。
SMO的伪代码如下:
创建一个alpha向量并将其初始化为零向量
当迭代次数小于最大迭代次数时(外循环)
对数据集中的每个特征向量(内循环):
如果该特征向量可以被优化:
随机选择另外一个特征向量
同时优化这两个向量
如果两个向量都不能被优化,退出内循环
如果所有向量都没有被优化,进行下一次迭代
python代码如下:
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):
#随机选择一个不等于i的j值
j=i
while (j==i):
j = int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L): #切一切alpha
#将alpha 限制在 L和 H之间
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
#SVM的一种实现,最小序列法(SMO)
#SMO 简化版
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 #用于记录alpha是否已经进行优化
for i in range(m): #对数据集遍历
fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b #预测的yi值(类别)
Ei = fXi - float(labelMat[i])#样本i的预测值与真值的误差
#检测样本是否违反KKT条件
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 #预测的yj值(类别)
Ej = fXj - float(labelMat[j]) #样本j的预测值与真值的误差
alphaIold = alphas[i].copy()#记录原始的alphai 的值
alphaJold = alphas[j].copy()#记录原始的alphaj 的值
#用于将aj调整到 0到 C之间
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:
continue
#eta 是 alphaj的最优修改量
eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
if eta >= 0: #作为简化的算法,这里省略了eta=0(很少发生)时求 alphaj的过程。
#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])#alphai 变化量和 alphaj 变化量 之和为0, 即保持 alphai+alphaj 不变
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):#成功的更新了一对alpha
iter_ += 1
else: iter_ = 0
#print ("iteration number: %d"% iter_)
return b,alphas
def calcWs(alphas,dataArr,classLabels):#计算权重系数
X = 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
def classfy(Xi, w, b): #做分类预测
y_predict = Xi*w + b
if y_predict>0:
return 1
else:
return -1
dataArr, labelArr = loadDataSet("testSet.txt")
#print(labelArr)
b, alphas = smoSimple(dataArr, labelArr, C=0.6, toler =0.001, maxIter=40)
w = calcWs(alphas, dataArr, labelArr)
print("偏置b: \n", b,"\n")
#print("alphas:\n", alphas,"\n")
print("权重矩阵w:\n", w,"\n") #预测值Y = Xi *w + b
print("支持向量:")
for i in range(len(alphas)):
if alphas[i]> 0: #打印支持向量
print(dataArr[i], labelArr[i])
print("第一个样本预测的分类:")
print(classfy(mat(dataArr)[0], w,b))
本文分享自 Python可视化编程机器学习OpenCV 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!