Machine Learning in Action -- Support Vector Machines

虽然SVM本身算法理论,水比较深,很难懂

但是基本原理却非常直观易懂,就是找到与训练集中支持向量有最大间隔的超平面

形式化的描述:

Machine Learning in Action -- Support Vector Machines

其中需要满足m个约束条件,m为数据集大小,即数据集中的每个数据点function margin都是>=1,因为之前假设所有支持向量,即离超平面最近的点,的function margin为1

对于这种有约束条件的最优化问题,用拉格朗日定理,于是得到如下的形式,

Machine Learning in Action -- Support Vector Machines

现在我们的目的就是求出最优化的m个拉格朗日算子,因为通过他们我们可以间接的算出w和b,从而得到最优超平面

考虑到某些离群值,加入松弛变量(slack variables),软间隔的版本为,

Machine Learning in Action -- Support Vector Machines

可以看到只是拉格朗日算子α的取值范围发生了变化,并且需要满足如下的KKT条件,

Machine Learning in Action -- Support Vector Machines

对于每个训练数据点都有一个对应的α
α对于绝大部分非支持向量而言,都是等于0的,这个由之前KKT条件可以推出
α等于C表示为离群值,个别function margin小于1的异常值
支持向量的α应该在0到C之间

所以这组KKT条件很容易理解,前面假设支持向量的function margin,即wx+b,为1

 

下面就是如何解这个问题?

之前大家都是用二次规划求解工具来求解,总之,这个计算量是非常大的

后来Platt提出SMO算法来优化这个求解过程,看过Andrew讲义就知道,这个算法是坐标上升算法的变种

下面就具体看看怎样实现SMO算法,关于算法的实现Andrew的讲义已经不够用,需要看Platt的论文

或参考这个blog,写的比较清楚,

支持向量机(五)SMO算法

 

SMO的简化版本

因为对于SMO算法,为了保持KKT条件,每次需要选取一对,即两个,α来进行优化
出于简单考虑,先随机的来选取α,后面再考虑启发式的来选取

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
b = 0; m,n = shape(dataMatrix) #初始化b为0
alphas = mat(zeros((m,1))) #初始化所有α为0
iter = 0
while (iter < maxIter):
alphaPairsChanged = 0
for i in range(m): #遍历所有的训练集
fXi = float(multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[i,:].T)) + b #1.计算wx+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(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]): #2.算出H和L,即α的边界值
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 - \ #3.计算新的α
dataMatrix[i,:]*dataMatrix[i,:].T - \
dataMatrix[j,:]*dataMatrix[j,:].T
if eta >= 0: print "eta>=0"; continue
alphas[j] -= labelMat[j]*(Ei - Ej)/eta #更新αj
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]) #根据αj的更新,反过来更新αi b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*\ #4.更新b
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): iter += 1 #记录α没有变化的次数,如果多次α没有变化,说明已经收敛
else: iter = 0
print "iteration number: %d" % iter return b,alphas

1. 计算wx+b

由于w是可以用α间接算出的,公式如下,

Machine Learning in Action -- Support Vector Machines

2. 计算α2的取值范围,即边界值,H和L,参考上面的blog

Machine Learning in Action -- Support Vector Machines

由约束条件,可以得到Machine Learning in Action -- Support Vector Machines

等式右边是常数,而y不是1,就是-1,所以α1和α2是上图的线性关系(y1=1,y2=-1的case)

图上画出线不同的两种边界值的case,可以比较直观的看出,

Machine Learning in Action -- Support Vector Machines

对于y1=1,y2=1的case同理差不多

3. 计算新的α

从讲义中,只知道算出下面式子的极值,就可以得到新的α

Machine Learning in Action -- Support Vector Machines

但没有给出具体解的过程,其实还是比较复杂的,参考上面的blog

得到的解是,

Machine Learning in Action -- Support Vector Machines

其中,

Machine Learning in Action -- Support Vector Machines K表示内积

看到这个公式就不难理解上面的代码的计算过程了

 

SMO的完整版本

关键的变化就是α的选择,采用启发式的来选择,从而大大优化算法速度

首先如何选取第一个α?

因为我们是需要通过最优化W来更新α的,所以对于边界值,即α=0或C的case,α是不会改变的,因为边界值的改变一定不会让分类效果变得更好

所以我们希望可以不断的优化支持向量的α,因为这个才是真正有效的优化,而不需要在边界α上浪费时间,这样可以大大提高算法效率

但问题是,所有α的初始值都是0,所以直接找非0的α是不可行的,所以这里的方法是,

先遍历所有α进行更新,这样会产生一些非0的α,然后对于这些α进行不断优化,直到收敛

然后再遍历所有α进行更新,也许会得到一些新的非0的α。。。。。

代码如下,

def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
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:
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:
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 #一次全遍历后,切换成选择遍历
elif (alphaPairsChanged == 0): entireSet = True #当前α不改变了,再全遍历 return oS.b,oS.alphas

上面的blog中,全遍历应该选取那些违反KKT条件的α进行优化,这样应该更合理一些,那些符合KTT条件的没必要做优化

有空可以把这个实现改一下,看看会不会更快一些

 

如何选择第二个α?

相对于之前的random选取,改进为,选择最大误差步长的,即Ei-Ej最大的

即当Machine Learning in Action -- Support Vector Machines为正时选择负的绝对值最大的Machine Learning in Action -- Support Vector Machines,反之,选择正值最大的Machine Learning in Action -- Support Vector Machines

这个直观上看,是有道理的,

因为另外一个α,需要根据当前α的改动,做相反的改动,以保证所有α的和不变

具体的代码参考原书,因为和简单版的相似

只是要缓存和更新所有的E值,并在选择第二个α时,做下计算选出Ei-Ej最大的

 

对非线性的复杂数据集使用kernels核函数

基础的SVM只能用于线性可分的数据集,如果不可分怎么办?

Machine Learning in Action -- Support Vector Machines

上面就是典型的例子,通常解决的方法,就是将低维数据转换为高维数据

比如上面的数据在二维看是线性不可分的,但是如果到了3维,可能蓝色的点和红色的点在高度上有些差异,这样就可以用一个平面,线性可分了

所以当数据非线性可分,当把数据的维度升高时,就有可能是线性可分的

对于SVM的对偶形式而言,只会计算内积形式,参考讲义

而我们如果把内积替换成核函数,就可以达成将低维数据到高维数据的转换,为什么?

因为所有的核函数,一定是等价于两个高维向量的内积,你不一定可以找出这个高维向量是什么,因为可能是无限维的向量,你也不需要care

所以用核函数替换当前低维向量的内积,就相当于用高维向量的内积替换低维向量的内积

下面看看具体的代码实现

这里实现的是最为流行的radial bias function的高斯版本,高斯核,对应于无限维的向量内积

Machine Learning in Action -- Support Vector Machines

对于m个数据点的训练集,我们要用高斯核函数,算出每对数据点间的核函数值,即m×m的结果矩阵

这样做训练或预测时,当需要计算内积时,直接查表即可

self.K = mat(zeros((self.m,self.m)))   #初始化m*m矩阵
for i in range(self.m): #对于一个数据点,算出和其他所有数据点之间的核函数值,填入表内
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

下面看下kernelTrans

def kernelTrans(X, A, kTup):
m,n = shape(X)
K = mat(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 = exp(K /(-1*kTup[1]**2))
else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
return K

其中kTup[0]表示核类型,这里给出两种核的实现

kTup[1]是计算核需要的参数,对于高斯核就是a

剩下使用时,就把所有需要计算内积的地方,都替换成查表即可

上一篇:一篇文章讲清楚android ImageView.ScaleType


下一篇:Codeforces Round #407 (Div. 2)A B C 水 暴力 最大子序列和