这个比较难懂,所以花费的时间比较长,我用了一个星期,希望看的人时间短一点。
先看https://blog.csdn.net/zouxy09/article/details/17291543,将他的思想结合下面的代码,看明白
还有这篇都结合着看 https://www.cnblogs.com/bentuwuying/p/6444249.html
对于H和L的理解那个我搞了半天,原理是这样
另外标签分类相等的一样的道理。
代码
python3版本代码(对于小白,在每次测试时候,可以打断点到测试那里,更容易理解,对于print过多,自行删除):
from random import random
import matplotlib.pyplot as plt
from numpy import *
#加载数据
def loadDataSet(fileName):
dataMat = []
label = []
fr = open(fileName)
for line in fr.readlines():
lineAr = line.strip().split('\t')
dataMat.append([float(lineAr[0]), float(lineAr[1])])
label.append(float(lineAr[2]))
return dataMat, label
#辅助函数1:随机一个不是第一个alpha的下标,为第二个alpha
def selectJrand(i, m):
"""
:param i:第一个alpha的下标
:param m: 所有alpha的数目
:return:随机一个不是第一个alpha下标为i的下标
"""
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
#辅助函数2:调整大于H,小于L的alpha,使aj落在L<aj<H这个区间上
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
#测试
dataArr, labelArr = loadDataSet('testSetSVM.txt')
print(labelArr)
#简化版的SMO,SMO表示序列最小优化(Sequential Minimal Optimization)。
#SMO工作原理:每次循环选取两个α,调整这两个值。
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
"""
:param dataMatIn: 数据集
:param classLabels: 类别标签
:param C: 常数C,C是离群点的权重,C越大表明离群点对目标函数影响越大,也就是越不希望看到离群点。
:param toler: 容错率->松弛变量
:param maxIter: 最大循环次数
:return:
"""
dataMatrix = mat(dataMatIn)
labelMat = mat(classLabels).T #转置为列向量
b = 0 #常数项b
m, n = shape(dataMatrix)
alphas = mat(zeros((m, 1)))
iter = 0 #迭代次数
while (iter < maxIter):
alphaPairsChanged = 0 #每次循环,都将他先设置为0,用于记录alpha是否优化
for i in range(m): #m为行
# 预测类别fXi = W.T * X + b,其中W =∑(alphas * label * X),下面表达式中后面dataMatrix[i, :].T是W里面的X,X表示传入的样本,下面的转置是矩阵计算问题
fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
Ei = fXi - float(labelMat[i]) #第一次误差Ei,如果误差很大就对该数据实例对应的alpha值进行优化
#检查是否有违反 KKT(加了松弛变量toler后) ,如果不满足,两种情况,以及优化
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
#这里选取另一个拉格朗日乘子aj,其他的拉格朗日乘子假设是定值AAA,由公式∑(alpha[n]*labelMat[n])=0,用它来表示ai*label[i]=0-AAA-aj*label[j]
j = selectJrand(i, m) #随机选择下一个alpha为aj
fXj = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b #对aj算预测
Ej = fXj - float(labelMat[j]) #再一次误差计算,为第二次误差Ej
alphaIold = alphas[i].copy() #记录他们以前的值
alphaJold = alphas[j].copy()
#算aj的取值范围
if (labelMat[i] != labelMat[j]): #不属于一个类别,求出aj的最大和最小取值范围
L = max(0, alphas[j] - alphas[i])
H = min(C, C + alphas[j] - alphas[i])
else: #属于同一类别,求出aj的最大和最小取值范围
L = max(0, alphas[j] + alphas[i] - C)
H = min(C, alphas[j] + alphas[i])
if L == H: #如果最底和最高相等,说明aj是个定值,可是我们要的aj是个变量,范围在L到H之间,不是我们想要的直接进行下次循环,本次循环结束
print("L==H")
continue
#eta是alpha[j]的最优修改量
eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
if eta >= 0: #如果eta>=0,直接进行下次循环,本次循环结束
print("eta>=0")
continue
#如果eta<0,更新alphas[j]
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
#与之前提到的一样aj[new]不是最终迭代后的值,需要再次进行约束:
alphas[j] = clipAlpha(alphas[j], H, L)
#检查alpha[j]是否有轻微的改变
if (abs(alphas[j] - alphaJold) < 0.00001):
print("j not moving enough")
continue
#这说明改变很大,更新alpha[i],修改量与j相同,但方向相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
#设置常数项b
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
#查找符合条件b1,b2,对b进行赋值
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))
# 如果没有优化,(就是没有进前面的for里面的if,也就是说,他们的点都在里面),查找下一个点看是否符合条件
if (alphaPairsChanged == 0):
iter += 1
else:
iter = 0
print("iteration number: %d" % iter)
return b, alphas
#测试
b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
print("由于aj的随机性每次运行都不一样:")
print("训练出来的b:", b)
print("大于0的alpha:", alphas[alphas > 0])
##############################################################################
#下面是画图部分,不想画就省略掉
#打印出支持向量(点)和W
print(shape(alphas), shape(mat(labelArr)), shape(dataArr))#看下形状,该转置转置等等
print("打印出的支持向量")
w = zeros((1, 2)) #w = ∑ (alphas * label * X),上面已经证明过,X为样本
for i in range(100):
if alphas[i] > 0.0: #拉格朗日乘子大于0的就是支持向量,等于0的为边缘线外的,对w没影响,没有小于0的
w += alphas[i].T*mat(labelArr[i]).T * dataArr[i] #因为这个dataArr[i]里面是全是x,y,所以算出来的w为
print(dataArr[i], labelArr[i])
print("打印出W:", w)
#进行画图
def plot(w, b):
"""
b+w[0]*x+w[1]*y=0
:param w:x和y前面的参数
:param b:就是截距
:return:
"""
dataMat, labelMat = loadDataSet('testSetSVM.txt') #加载数据集,返回数据集合分类
dataArr = array(dataMat) #将数据集数组化
dataShape = shape(dataArr) #data的形状
n = dataShape[0] #形状的第一个是行数(即数据集的个数),第二个是列数(即数据集的特征)
xcord1 = [] #分类为1的点
ycord1 = []
xcord0 = [] #分类为-1的点
ycord0 = []
for i in range(n):
if int(labelMat[i]) == 1: #如果分类为1,添加到1分类的点集,否者返回到-1分类的点集
xcord1.append(dataArr[i, 0])
ycord1.append(dataArr[i, 1])
else:
xcord0.append(dataArr[i, 0])
ycord0.append(dataArr[i, 1])
fig = plt.figure(figsize=(6, 6)) #figure()函数创建一个新的图
ax = fig.add_subplot(111) #add_subplot()函数在一张figure里面生成多张子图参数111,表示1行1列第1个位置
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s') #散点图
ax.scatter(xcord0, ycord0, s=30, c='green')
#对支持向量画圈
for i in range(100):
if alphas[i] > 0.0:
ax.scatter(dataArr[i, 0], dataArr[i, 1], color='', marker='o', edgecolors='b', s=200)
# 直线 x 的范围
x = arange(-2.0, 10.0, 0.1) # (start, end, step)可以调step执行起来,看看图
# 画出直线,b+w[0]*x+w[1]*y=0
# 实际他是一个 y=ax+b, b常量
y = (-b[0]-w[0]*x)/w[1] #这个就是我们求出来的模型函数
ax.plot(x, y) #画出直线
plt.xlabel('x1') #X轴标记为x1
plt.ylabel('x0') #Y轴标记为x0
plt.show()
#这个传值进去的必须是数组.不是数组的话,画线的时候y=...就不会被broadcast,
# w为数组,b为矩阵,所以将b用getA()转为数组
# w的转置,因为原本shape(w)是(1,2)的,他转化为数组不能分为w[0],w[1],画线的时候会报错,ValueError: operands could not be broadcast together with shapes (2,) (120,)
plot(w.T, b.getA())
testSetSVM.txt 测试集
3.542485 1.977398 -1
3.018896 2.556416 -1
7.551510 -1.580030 1
2.114999 -0.004466 -1
8.127113 1.274372 1
7.108772 -0.986906 1
8.610639 2.046708 1
2.326297 0.265213 -1
3.634009 1.730537 -1
0.341367 -0.894998 -1
3.125951 0.293251 -1
2.123252 -0.783563 -1
0.887835 -2.797792 -1
7.139979 -2.329896 1
1.696414 -1.212496 -1
8.117032 0.623493 1
8.497162 -0.266649 1
4.658191 3.507396 -1
8.197181 1.545132 1
1.208047 0.213100 -1
1.928486 -0.321870 -1
2.175808 -0.014527 -1
7.886608 0.461755 1
3.223038 -0.552392 -1
3.628502 2.190585 -1
7.407860 -0.121961 1
7.286357 0.251077 1
2.301095 -0.533988 -1
-0.232542 -0.547690 -1
3.457096 -0.082216 -1
3.023938 -0.057392 -1
8.015003 0.885325 1
8.991748 0.923154 1
7.916831 -1.781735 1
7.616862 -0.217958 1
2.450939 0.744967 -1
7.270337 -2.507834 1
1.749721 -0.961902 -1
1.803111 -0.176349 -1
8.804461 3.044301 1
1.231257 -0.568573 -1
2.074915 1.410550 -1
-0.743036 -1.736103 -1
3.536555 3.964960 -1
8.410143 0.025606 1
7.382988 -0.478764 1
6.960661 -0.245353 1
8.234460 0.701868 1
8.168618 -0.903835 1
1.534187 -0.622492 -1
9.229518 2.066088 1
7.886242 0.191813 1
2.893743 -1.643468 -1
1.870457 -1.040420 -1
5.286862 -2.358286 1
6.080573 0.418886 1
2.544314 1.714165 -1
6.016004 -3.753712 1
0.926310 -0.564359 -1
0.870296 -0.109952 -1
2.369345 1.375695 -1
1.363782 -0.254082 -1
7.279460 -0.189572 1
1.896005 0.515080 -1
8.102154 -0.603875 1
2.529893 0.662657 -1
1.963874 -0.365233 -1
8.132048 0.785914 1
8.245938 0.372366 1
6.543888 0.433164 1
-0.236713 -5.766721 -1
8.112593 0.295839 1
9.803425 1.495167 1
1.497407 -0.552916 -1
1.336267 -1.632889 -1
9.205805 -0.586480 1
1.966279 -1.840439 -1
8.398012 1.584918 1
7.239953 -1.764292 1
7.556201 0.241185 1
9.015509 0.345019 1
8.266085 -0.230977 1
8.545620 2.788799 1
9.295969 1.346332 1
2.404234 0.570278 -1
2.037772 0.021919 -1
1.727631 -0.453143 -1
1.979395 -0.050773 -1
8.092288 -1.372433 1
1.667645 0.239204 -1
9.854303 1.365116 1
7.921057 -1.327587 1
8.500757 1.492372 1
1.339746 -0.291183 -1
3.107511 0.758367 -1
2.609525 0.902979 -1
3.263585 1.367898 -1
2.912122 -0.202359 -1
1.731786 0.589096 -1
2.387003 1.573131 -1