https://blog.csdn.net/sinat_20177327/article/details/79729551
https://blog.csdn.net/c406495762/article/details/78072313
https://blog.csdn.net/willbkimps/article/details/54697698
数据源地址:https://github.com/pbharrin/machinelearninginaction/blob/master/Ch06/testSet.txt
目录
1、拉格朗日函数和拉格朗日乘子
要找到函数 z=f(x,y) 在附加条件 φ(x,y)=0 ,下可能的极值点,可以先引进拉格朗日函数:
其中 λ 称为拉格朗日乘子。求其对 x 与 y 的一阶偏导数,并使之为零,然后再于附加条件方程联立:
同理可以扩展至多元函数:
2、超平面:这里说的超平面实际为仿射超平面
仿射超平面是仿射空间中的代数1的仿射子空间。 在笛卡尔坐标中,可以用以下形式的单一线性方程来描述这样的超平面(至少有一个ai不是0):
在真实仿射空间的情况下,换句话说,当坐标为实数时,该仿射空间将空间分成两个半空间,它们是超平面的补码的连接分量,由不等式给出:
和
作为一个例子,一点是一维空间中的超平面,一条线是二维空间中的超平面,一个平面是三维空间中的超平面。 三维空间中的一行不是超平面,并且不将空间分成两部分(这样一条线的补码连接起来)。
欧几里德空间的任何超平面具有两个单位法向量。
仿射超平面用于定义许多机器学习算法中的决策边界,例如线性组合(倾斜)决策树和感知器。
3、凸函数:
首先,以上的支持向量机是凸函数。理解凸函数,我们还要先明确另一个概念,凸集。在凸几何中,凸集(convex set)是在)凸组合下闭合的放射空间的子集。看一幅图可能更容易理解:
左右量图都是一个集合。如果集合中任意2个元素连线上的点也在集合中,那么这个集合就是凸集。显然,上图中的左图是一个凸集,上图中的右图是一个非凸集。
凸函数的定义也是如此,其几何意义表示为函数任意两点连线上的值大于对应自变量处的函数值。若这里凸集C即某个区间L,那么,设函数f为定义在区间L上的函数,若对L上的任意两点x1,x2和任意的实数λ,λ属于(0,1),总有:
则函数f称为L上的凸函数,当且仅当其上镜图(在函数图像上方的点集)为一个凸集。再看一幅图,也许更容易理解:
像上图这样的函数,它整体就是一个非凸函数,我们无法获得全局最优解的,只能获得局部最优解。比如红框内的部分,如果单独拿出来,它就是一个凸函数
4、线性可分支持向量机:
首先用向量形式表达超平面方程:
其中w为法向量,决定了超平面的方向,b为位移量,可以通过三维坐标系的点到平面的距离来理解点到超平面的距离,因为点到超平面的距离没有物理概念,只有数学概念:
点到平面的距离: 点(x0,y0,z0)到了平面Ax+By+Cz+D=0的距离 为:d=|Ax0+By0+Cz0+D|/√(A^2+B^2+C^2)
假设超平面能将训练样本正确地分类,即对于训练样本,满足以下公式:
公式称为最大间隔假设,yi=+1表示样本为正样本,yi=−1 表示样本为负样本,式子前面选择大于等于+1,小于等于-1只是为了计算方便,原则上可以是任意常数,但无论是多少,都可以通过对 w 的变换使其为 +1 和 -1 ,此时将公式左右都乘以 yi,得到如下:
实际上等价于(训练集中的所有样本都应满足的公式):
如下图所示,距离超平面最近的这几个样本点满足 ,它们被称为“支持向量”。虚线称为边界,两条虚线间的距离称为间隔(margin):
下面我们开始计算间隔,其实间隔就等于两个异类支持向量的差在 w 上的投影,即:
其中 x+ 和 x−分别表示两个正负支持向量,因为 x+ 和 x−满足 ,即:
推出:
代入w 上的投影公式中可以得到:
至此,我们求得了间隔,SVM的思想是使得间隔最大化,也就是:
显然,最大化 2/||w||相当于最小化 ||w||,为了计算方便,将公式转化成如下:
以上公式即为支持向量机的基本型。
5、线性可分支持向量机求极值的拉格朗日函数转化:
我们通过拉格朗日乘子来求解:
分别对 w 和 b求偏导:
令其分别为0,可以得到:
将公式代入拉格朗日函数中,可得:
此时,原问题就转化为以下仅关于 α 的问题:
解出 α 之后,根据可以求得 w , 进而求得 b,可以得到模型:
上述过程的KKT条件为:
我们分析一下,对于任意的训练样本 (xi,yi),
- 若 αi=0,则其不会在模型中的求和项中出现,也就是说,它不影响模型的训练;
- 若 αi>0,则 yif(xi)−1=0,也就是 yif(xi)=1,即该样本一定在边界上,是一个支持向量。
这里显示出了支持向量机的重要特征:当训练完成后,大部分样本都不需要保留,最终模型只与支持向量有关。
6、松弛变量:
至此,一切都很完美,但是这里有个假设:数据必须100%线性可分。目前为止,我们知道 几乎所有数据都不那么“干净”。这时我们就可以通过引入所谓松弛变量(slack variable),来允 许有些数据点可以处于分隔面的错误一侧。这样我们的优化目标就能保持仍然不变,但是此时新 的约束条件则变为:
这里的常数C用于控制“大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的 权重。在优化算法的实现代码中,常数C是一个参数,因此我们就可以通过调节该参数得到不同 的结果。一旦求出了所有的alpha,那么分隔超平面就可以通过这些alpha来表达。这一结论十分 直接,SVM中的主要工作就是求解这些alpha。
7、对偶问题:
对于一个问题A:找到函数f(x)的最小值的x,
可以找到对应的另外一个问题B:找到函数u(x)的最大值的x,
并且,B中的x就是A中的x,
那么,A和B就可以理解为对偶问题。
8、简化版SMO算法
现在,我们的问题转变成求解式5中提及的拉格朗日函数。这是一个二次规划问题,可以使用通用的二次规划算法来求解。这里可以使用更高效的算法,SMO算法。
SMO的基本思路是先固定之外的所有参数,然后求上的极值。由于存在约束,若固定之外的其他参数,则可由其他变量导出。于是,SMO每次选择两个变量和,并固定其他参数。这样,在参数初始化后,SMO不断执行如下两个步骤直至收敛:
- 选取一对需要更新的变量和
- 固定和以外的参数,求解拉格朗日函数获得更新后的和
在这里我们假设正在优化的两个拉格朗日乘子对应的样本正好是第一个和第二个,两个拉格朗日乘子分别为和,他们满足的约束可以写为:
其中,
c是常量,这样就固定了,i≥3。我们可以用表示,这样问题就转化成为关于 的一元二次函数优化问题,然后对求极值。确定下来了,也就确定了。现在,我们先求的上下限。
和 的值为+1或者-1。和 的约束条件式可以用下图表示,横轴为,纵轴为:
从图中可以看出的约束条件:
当≠ 时,L=max(0, -),H=min(C, C+- )
当= 时,L=max(0, +),H=min(C, +-C)
因此有:
将和代入拉格朗日函数,求最大值:
下图为K(, )的每一项,用表示K(, )
所以:
和④都是常量,因此记ψconstant = + ④
将式f(x)变形得:
记:
因此:
又因为:
两边同时乘以,得:
记:
那么:
代入ψ(α ),因此有:
对求偏导有:
如果二阶导数小于0,当一阶导数为0时有极大值。令上式等于0,有:
将h和代入继续推导:
即有:
记:
有:
b的更新:
Python实现:
# coding:UTF-8
from numpy import *
from time import sleep
from matplotlib.patches import Circle
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 #we want to select any J not equal to 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): #五个输入参数对应表示数据集,类别标签,常数C,容错率,退出前最大循环次数
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose() #将输入数据转化为matrix的形式,将类别标签进行转置
b = 0; m,n = shape(dataMatrix) #将常数b进行初始化,m,n分别表示数据集的行和列
alphas = mat(zeros((m,1))) #将alpha初始化为m行1列的全零矩阵
iter = 0 #该变量用于存储是在没有任何alpha改变的情况下遍历数据集的次数
while (iter < maxIter): #当遍历数据集的次数超过设定的最大次数后退出循环
alphaPairsChanged = 0 #该参数用于记录alpha是否已进行优化,每次循环先将其设为0
for i in range(m): #将数据集中的每一行进行测试
#multiply是numpy中的运算形式,将alpha与标签相乘的结果的转置与对应行的数据相乘再加b
#fxi表示将该点带入后进行计算得到的预测的类别,若Ei=0,则该点就在回归线上
fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
#label[i]*Ei代表误差,如果误差很大可以根据对应的alpha值进行优化,同时检查alpha值使其不能等于0或C
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m) #随机选取另外一个数作为j,也就是选取另外一行数据
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j]) #采用与上述相同的方式来计算j行的预测误差
alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy(); #暂存变量,便于后面对其进行比对
#调整H与L的数值,便于用来控制alpha的范围
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 #continue在python中代表本次循环结束直接运行下一次循环
#eta代表的是alphas[j]的最优修改量,如果eta大于等于0则需要退出for循环的当前迭代过程
eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
if eta >= 0: continue
#给alphas[j]赋新值,同时控制alphas[j]的范围在H与L之间
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
alphas[j] = clipAlpha(alphas[j],H,L)
#然后检查alphas[j]是否有轻微改变,如改变很小,就进行下一次循环
if (abs(alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); continue
#将alpha[i]与alpha[j]同时进行改变,改变的数值保持一致,但是改变的方向相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
#在完成alpha[i]与alpha[j]的优化之后,给这两个alpha设置一个常数项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
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 #该参数用于记录alpha是否已进行优化
print ("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) #iter表示遍历的次数
if (alphaPairsChanged == 0): 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 plotBestFit(dataMatIn,labelMatIn,alphas,ws_float0,ws_float1,b_float):
import matplotlib.pyplot as plt
dataMat,labelMat=loadDataSet('F://电子书/testSetSVN.txt')
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
for i in range(n):
if int(labelMat[i])== 1:
xcord1.append(dataArr[i,0]); ycord1.append(dataArr[i,1])
else:
xcord2.append(dataArr[i,0]); ycord2.append(dataArr[i,1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
plotx = arange(0.0 , 10.0 , 0.1)
ploty = -(ws_float0*plotx+ b_float)/ws_float1
ax.plot(plotx,ploty)
for i in range(100):
if alphas[i]>0.0:
circle = Circle((dataArr[i,0], dataArr[i,1]), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
ax.add_patch(circle)
plt.xlabel('X1'); plt.ylabel('X2');
plt.show()
if __name__ == '__main__':
dataMat,labelMat = loadDataSet('F://电子书/testSetSVN.txt')
b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
print ('b:',b)
print ('alphas[alphas>0]:',alphas[alphas>0])
print (shape(alphas[alphas>0]))
for i in range(100):
if alphas[i]>0.0:
print (dataMat[i],labelMat[i])
ws = calcWs(alphas,dataMat,labelMat)
b_float = float(b[0])
ws_float0 = float(ws[0])
ws_float1 = float(ws[1])
plotBestFit(dataMat,labelMat,alphas,ws_float0,ws_float1,b_float);
9、完整的Platt SMO算法
主要加速在于smoP方法中,不是每次迭代都去便利所有的数据集,一开始的for循环在数据集上遍历任意可能 的alpha 。我们通过调用innerL()来选择第二个alpha,并在可能时对其进行优化处理。如果有 任意一对alpha值发生改变,那么会返回 1。第二个for循环遍历所有的非边界alpha值,也就是不 在边界0或C上的值。
Python实现:
# coding:UTF-8
from numpy import *
from time import sleep
from matplotlib.patches import Circle
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 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 selectJ(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 selectJrand(i,m):
j=i #we want to select any J not equal to 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 updateEk(oS, k):#after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
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 #changed for kernel
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
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 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
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 plotBestFit(dataMatIn,labelMatIn,alphas,ws_float0,ws_float1,b_float):
import matplotlib.pyplot as plt
dataMat,labelMat=loadDataSet('F://电子书/testSetSVN.txt')
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
for i in range(n):
if int(labelMat[i])== 1:
xcord1.append(dataArr[i,0]); ycord1.append(dataArr[i,1])
else:
xcord2.append(dataArr[i,0]); ycord2.append(dataArr[i,1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
plotx = arange(0.0 , 10.0 , 0.1)
ploty = -(ws_float0*plotx+ b_float)/ws_float1
ax.plot(plotx,ploty)
for i in range(100):
if alphas[i]>0.0:
circle = Circle((dataArr[i,0], dataArr[i,1]), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
ax.add_patch(circle)
plt.xlabel('X1'); plt.ylabel('X2');
plt.show()
if __name__ == '__main__':
dataMat,labelMat = loadDataSet('F://电子书/testSetSVN.txt')
b,alphas = smoP(dataMat, labelMat, 0.6, 0.001, 40)
print ('b:',b)
print ('alphas[alphas>0]:',alphas[alphas>0])
print (shape(alphas[alphas>0]))
for i in range(100):
if alphas[i]>0.0:
print (dataMat[i],labelMat[i])
ws = calcWs(alphas,dataMat,labelMat)
b_float = float(b[0])
ws_float0 = float(ws[0])
ws_float1 = float(ws[1])
plotBestFit(dataMat,labelMat,alphas,ws_float0,ws_float1,b_float);
10、非线性支持向量机和核函数:
对于非线性问题,线性可分支持向量机并不能有效解决,要使用非线性模型才能很好地分类。先看一个例子,如下图,很显然使用直线并不能将两类样本分开,但是可以使用一条椭圆曲线(非线性模型)将它们分开。非线性问题往往不好求解,所以希望能用解线性分类问题的方法求解,因此可以采用非线性变换,将非线性问题变换成线性问题。
对于这样的问题,可以将训练样本从原始空间映射到一个更高维的空间,使得样本在这个空间中线性可分,如果原始空间维数是有限的,即属性是有限的,那么一定存在一个高维特征空间是样本可分。令ϕ(x)表示将 x 映射后的特征向量,于是在特征空间中,划分超平面所对应的的模型可表示为:
于是有最小化函数:
其对偶问题为:
若要对公式(16)求解,会涉及到计算 ,这是样本和 映射到特征空间之后的内积,由于特征空间的维数可能很高,甚至是无穷维,因此直接计算 通常是困难的,于是想到这样一个函数:
即 和 在特征空间中的内积等于他们在原始样本空间中通过函数计算的函数值,于是公式(16)写成如下:
求解后得到:
这里的函数 就是核函数,在实际应用中,通常人们会从一些常用的核函数里选择(根据样本数据的不同,选择不同的参数,实际上就得到了不同的核函数),下面给出常用的核函数:
- 线性核:
- 多项式核(d是多项式的次数,d=1是退化为线性核):
- 高斯核(σ>0):
- 拉普拉斯核(σ>0):
- sigmiod核(β>0,θ>0):
此外,核函数也可以通过组合得到,在此不再赘述。