码迷,mamicode.com
首页 > 其他好文 > 详细

Machine Learning in Action -- Support Vector Machines

时间:2014-07-24 21:45:22      阅读:332      评论:0      收藏:0      [点我收藏+]

标签:style   blog   http   color   使用   os   数据   io   

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

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

形式化的描述:

bubuko.com,布布扣

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

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

bubuko.com,布布扣

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

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

bubuko.com,布布扣

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

bubuko.com,布布扣

对于每个训练数据点都有一个对应的α
α对于绝大部分非支持向量而言,都是等于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是可以用α间接算出的,公式如下,

bubuko.com,布布扣

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

bubuko.com,布布扣

由约束条件,可以得到
bubuko.com,布布扣

等式右边是常数,而y不是1,就是-1,所以α1和α2是上图的线性关系(y1=1,y2=-1的case)
图上画出线不同的两种边界值的case,可以比较直观的看出,
bubuko.com,布布扣

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

3. 计算新的α

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

bubuko.com,布布扣

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

得到的解是,

bubuko.com,布布扣

其中,

bubuko.com,布布扣 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最大的

即当bubuko.com,布布扣为正时选择负的绝对值最大的bubuko.com,布布扣,反之,选择正值最大的bubuko.com,布布扣

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

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

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

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

 

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

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

bubuko.com,布布扣

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

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

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

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

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

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

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

下面看看具体的代码实现

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

bubuko.com,布布扣

对于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

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

Machine Learning in Action -- Support Vector Machines,布布扣,bubuko.com

Machine Learning in Action -- Support Vector Machines

标签:style   blog   http   color   使用   os   数据   io   

原文地址:http://www.cnblogs.com/fxjwind/p/3865990.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!