支持向量机-（1）

怼烎@ 2022-06-02 07:18 196阅读 0赞

要开始SVM和SMO之前需要掌握以下几个概念：

1超平面

将数据分割成两部分，一部分\-1，一部分\+1，最优的分割方式这个公式被称为超平面，n维的超平面是n-1维的，所以2维的超平面是线，3维的超平面是面，以此类推。

超平面表现方式为：f(x)=WX+b，求超平面的过程也就是求W和b的过程，其中W称为法向量，b称为截距。

![20180107223139326][]

2支持向量

距离超平面距离最近的两边的点是支持向量，两边的支持向量距离超平面相等且最近。

![20180107223208930][]

3松弛变量

并不是所有的数据都是线性可分的，为了使得非线性的情况也可以被线性处理，引入了松弛变量的概念。引入松弛变量后只是将线性的概率大大增加，并不能保证一定可以线性化。

4函数间隔

在超平面 WX+b=0 确定的情况下， | Wx+b | 标识点距离超平面的远近， Wx+b 结果与 y 的符号是否一致代表分类是否准确，所以将 y\*( Wx+b ) 定义为函数间隔，代表分类的确信度和准确读。简写为γ=y\*f(x)

5几何间隔

在选择超平面时只通过函数间隔是不足够的，因为例如NW\*X+Nb这种函数间隔等比例增大后，超平面还是不变的，所以需要引入几何间隔的概念。引入||w||=(W\*\*2+b\*\*2)\*\*0.5这个缩放因子，几何间隔γ=y\*f(x)/||w||，所有的点γp>=γ，其中相等的点就是支持向量。如果我们认为支持向量的γ=1，那么取几何间隔最大化的问题就演变成max(1/||w||)的问题，前提条件是y\*f(x)>=1

6对偶算法

原始问题比较难求解时求解另外一个问题，例如我们无法求解如何使|X|最大时，我们可以求解如何使1/|x|最小，这两个问题是对偶的

7拉格朗日乘值

通过拉格朗日函数将约束条件融合到目标函数里去，融合后的结果是条件消失，函数增加一个α参数，例如5中几何间隔最大化的问题就是典型的符合拉格朗日函数的例子，融合后y\*f(x)>=1这个条件消失。

8矩阵的转置：

![20180107223352453][]

推论的过程：

第一步，引入拉格朗日：

通过前面几何间隔最大化已知我们要求解的是求 max(1/||w||) 在 y\*f(x)>=1 条件下的最大化，通过对偶，转化成求![20180107223427728][]，条件不变，引入拉格朗日乘值以后条件消失，函数变成：

![20180107223459441][]

由于原来的前提条件是y\*f(x)>=1，所以可以推断的出括号内>=0，所以函数里α>=0

(PS:不是说拉格朗日是消除条件的么，结果是又带进来新的条件，不过看在这个新的条件比较简单的份子上，我们就原谅他了)

因为||w||也是由w和b组成的，那么现在的问题就变成了求w、b、α三者组合的问题

![20180107223549896][]

先求对w，b的极小，再求对α的极大。

第二步，求偏导做转化：

先假设α固定，对w和b求偏导

![20180107223632192][]

带入第一步的函数后演算：

![20180107223658499][]

最终演化成：

![20180107223727940][]

那么这时候惊喜的发现整个函数中所有y和x都是我们训练集里的已知数，只剩下α这一个变量！。

SMO介绍：

![20180107223803954][]

是否有更快捷的解法？

关于SMO的推论网上有个很好的教材[https://www.jianshu.com/p/55458caf0814][https_www.jianshu.com_p_55458caf0814]，需要慢慢加强理解和推演，这里先介绍如何直接使用：

1先给定假设的α，然后算出2个预测结果f(x)1、f(x)2，并根据结果集算出误差E1=y1-f(x)1, E2=y2-f(x)2,

![20180107223846840][]

2根据y1、y2是否同号，给定α的上限和下限

![20180107223916753][]

3计算α2的新值

![20180107223946440][]

其中

![20180107224015345][]

4根据2中L、H的限制，修订α2

![20180107224059081][]

5更新α1

![20180107224131952][]

6更新b

![20180107224202017][]

![20180107224227491][]

7不断循环优化α和b，收敛后根据α求w

![20180107224259751][]

8算出w和b之后α也就没用了，根据w、b对数据做预测，预测结果大于0则为1，小于0则为\-1

f(x)=wx+b

本篇先介绍机器学习实践中可以直接被线性处理的这个案例，训练集在git上(https://github.com/yejingtao/forblog/blob/master/MachineLearning/trainingSet/svm/testSet.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

基础的辅助函数：

#很简单，加载训练集，返回特性集和结果集
    def loadDataSet (fileName) :
        dataMat = []
        labelMat = []
        fr = open(fileName)
        for line in fr.readlines() :
            lineArray = line.strip().split('\t')
            dataMat.append([float(lineArray[0]),float(lineArray[1])])
            labelMat.append(float(lineArray[2]))
        return dataMat, labelMat
    
    #返回一个0-m之间非i的随机数
    def selectJrand (i, m) :
        j = i
        while j==i :
            j = int(random.uniform(0,m))
        return j
    
    def clipAlpha (aj,H,L) :
        if aj>H : return H
        if aj<L : return L
        return aj

下面是核心的处理函数：

#C是松弛度
    #toler是容错率
    #maxInter迭代次数
    def smoSimple(dataMatIn,labelMatIn, C, toler, maxInter) :
        dataMat = mat(dataMatIn)
        labelMat = mat(labelMatIn).transpose()
        #初始化b=0
        b=0
        m,n = shape(dataMat)
        #初始化固定的alpha
        alphas = mat(zeros((m,1)))
        iter=0
        while iter<maxInter :
            isAlphaChanged = False
            for i in range(m) :
                #基于目前α给出的预测结果
                fXi = float(multiply(alphas,labelMat).T * (dataMat * dataMat[i,:].T)) + b
                #预测的误差
                Ei = fXi - float(labelMat[i])
                #如果误差很大，说明α需要优化，进入优化逻辑
                if ((labelMat[i]*Ei<-toler) and (alphas[i]<C)) or \
                        ((labelMat[i] * Ei > toler) and (alphas[i] > 0)) :
                    #α由于存在yα积分为0的约束，所以必须成对的进行优化，随机选一个非i的j
                    j = selectJrand(i,m)
                    #计算使用αj的预测结果
                    fXj = float(multiply(alphas,labelMat).T * (dataMat * dataMat[j,:].T)) + b
                    #计算αj预测的误差
                    Ej = fXj - float(labelMat[j])
                    #传值方式备份一次需要优化的αi和αj
                    alphaOldi = alphas[i].copy()
                    alphaOldj = alphas[j].copy()
                    #如果y不同号，求L和H
                    if labelMat[i] != labelMat[j] :
                        L = max(0,alphas[j]-alphas[i])
                        H = min(C, C+alphas[j]-alphas[i])
                    # 如果y同号，求L和H
                    else :
                        L = max(0,alphas[j]+alphas[i]-C)
                        H = min(C, alphas[j]+alphas[i])
                    if L==H :
                        #print('L=H')
                        continue
                    #求j的优化系数eta
                    eta = 2.0*dataMat[i,:]*dataMat[j,:].T - \
                          dataMat[i, :] * dataMat[i, :].T - \
                          dataMat[j, :] * dataMat[j, :].T
                    #为啥大于零就退出了？
                    if eta >=0 :
                        print('eta>0')
                        continue
                    #优化αj
                    alphas[j] -= labelMat[j] * (Ei-Ej)/eta
                    #限制α范围
                    alphas[j] = clipAlpha(alphas[j],H,L)
                    #优化步长不够跳过优化，从某种程度来说0.0001可以代表选取超平面的精度
                    if (abs(alphas[j]-alphaOldj)<0.0001) :
                        #print('j not moving enough')
                        continue
                    #更新αi
                    alphas[i] += labelMat[j] * labelMat[i] * (alphaOldj-alphas[j])
                    #以下优化b
                    b1 = b - Ei - labelMat[i]*(alphas[i]-alphaOldi)*dataMat[i, :] * dataMat[i, :].T \
                        -labelMat[j]*(alphas[j]-alphaOldj)*dataMat[i, :] * dataMat[j, :].T
                    b2 = b - Ej - labelMat[i] * (alphas[i] - alphaOldi) * dataMat[i, :] * dataMat[j, :].T \
                         - labelMat[j] * (alphas[j] - alphaOldj) * dataMat[j, :] * dataMat[j, :].T
                    if (0<alphas[i]) and (alphas[i]<C) : b = b1
                    elif (0<alphas[j]) and (alphas[j]<C) : b = b2
                    else : b = (b1+b1)/2.0
                    isAlphaChanged = True
            if isAlphaChanged==False : iter+=1
            else : iter=0
        return b,alphas

为什么α默认矩阵是0？

因为只有这样才能保证:

![20180108132852226][]  
上面不断收敛最终求得b和alphas，下面再根据alphas求得w

#根据alpha计算出w
    def calcWs(alpha, dataArr, labelArr) :
        dataMat = mat(dataArr)
        labelMat = mat(labelArr).transpose()
        m,n = shape(dataMat)
        w = zeros((n,1))
        for i in range(m) :
            w += multiply(alpha[i]*labelMat[i],dataMat[i,:].T)
        return w

虽然上述for循环遍历了数据集中的所有数据，但是最终起作用的只有支持向量。

现在w和b都全了，下面可以进行预测了：

#利用w和b预测结果
    def validate(dataMat, labelArr, w, b):
        m,n = shape(dataMat)
        for i in range(m) :
            y = dataMat[i,:]*mat(w)+b
            if y>0 :
                resultY=1
            else :
                resultY = -1
            if labelArr[i]==resultY :
                print('%d is right, result is %f' %(i,y))
            else :
                print('%d is error, result is %f' %(i,y))

测试代码：

dataMat, labelMat = loadDataSet('C:\\2017\\提高\\机器学习\\训练样本\\svm\\testSet.txt')
    b,alphas = smoSimple(dataMat,labelMat,0.6,0.001,40)
    print(b)
    print(alphas[alphas>0])
    #[[-3.82650687]]
    #[[ 0.12688066  0.24039407  0.36727473]]
    w = calcWs(alphas, dataMat, labelMat)
    print(w)

对上面代码进行改造后：

class optStruts :
        def __init__(self,dataMatIn,labelMatIn, C, toler):
            self.X = dataMatIn
            self.Y = labelMatIn
            self.C = C
            self.toler = toler
            self.m = shape(dataMatIn)[0]
            self.alphas = mat(zeros((self.m,1)))
            self.b = 0
            #误差缓存，第一列为是否有效（0无效，1有效），第二列为缓存的值
            self.eCache = mat(zeros((self.m,2)))
    
    #求k列的误差
    def calcEk(os, k):
        fXk  = float(multiply(os.alphas,os.Y).T * (os.X * os.X[k,:].T)) + os.b
        #误差
        ek = fXk - os.Y[k]
        return ek
    
    #用于选择第二个α，也就是αj
    #保证选择的αj可以优化最大的步长，参考前面的公式labelMat[j] * (Ei-Ej)/eta
    def selectJ(i, os, Ei):
        maxK = -1
        maxDeleteE = 0
        Ej = 0
        #1标识设置缓存为有效
        os.eCache[i] = [1,Ei]
        #构造非零E值所对应的alpha值构造的列表
        validateEcacheList = nonzero(os.eCache[:,0].A)[0]
        if len(validateEcacheList) >1 :
            for k in validateEcacheList :
                #配对优化，不能优化自己
                if k==i :continue
                #挨个算误差
                Ek = calcEk(os, k)
                deltaE = abs(Ei-Ek)
                if deltaE>maxDeleteE :
                    maxK = k
                    maxDeleteE = deltaE
                    Ej = Ek
            #返回具有最大误差的α坐标和误差Ek
            return maxK, Ej
        else :
            #都为零，说明是起步阶段，随便选个非i的运算即可。
            j = selectJrand(i,os.m)
            Ej = calcEk(os, j)
        return j, Ej
    
    #算k的误差，并设置为有效
    def updateEk(os , k) :
        Ek = calcEk(os, k)
        os.eCache[k] = [1, Ek]
    
    def innerL (i, os) :
        #获取i的误差
        Ei = calcEk(os,i)
        if ((Ei*os.Y[i]<-os.toler) and (os.alphas[i]<os.C)) or \
                ((Ei*os.Y[i]>os.toler) and (os.alphas[i]>0)) :
            j, Ej = selectJ(i, os, Ei)
            alphas_i_old = os.alphas[i].copy()
            alphas_j_old = os.alphas[j].copy()
            if os.Y[i] != os.Y[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 H==L : 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
            if eta>=0 : return 0
            os.alphas[j] -= os.Y[j]*(Ei-Ej)/eta
            os.alphas[j] = clipAlpha(os.alphas[j],H,L)
            updateEk(os,j)
            if (abs(os.alphas[j]-alphas_j_old)<0.0001) : return 0
            os.alphas[i] += os.Y[i]*os.Y[j] * (alphas_j_old - os.alphas[j])
            updateEk(os, i)
            b1 = os.b - Ei - os.Y[i] * (os.alphas[i] - alphas_i_old) * os.X[i, :] * os.X[i, :].T \
                 - os.Y[j] * (os.alphas[j] - alphas_j_old) * os.X[i, :] * os.X[j, :].T
            b2 = os.b - Ej - os.Y[i] * (os.alphas[i] - alphas_i_old) * os.X[i, :] * os.X[j, :].T \
                 - os.Y[j] * (os.alphas[j] - alphas_j_old) * os.X[j, :] * os.X[j, :].T
            if (0 < os.alphas[i]) and (os.alphas[i] < os.C):
                os.b = b1
            elif (0 < os.alphas[j]) and (os.alphas[j] < os.C):
                os.b = b2
            else:
                os.b = (b1 + b1) / 2.0
            return 1
        else :
            return 0
    
    def smop(dataMatIn,labelMatIn, C, toler, maxInter, kTup=('lin',0)) :
        #构造os辅助对象
        os = optStruts(mat(dataMatIn),mat(labelMatIn).transpose(), C, toler)
        iter = 0
        entireSet = True
        alphaPairsChanged = 0
        while (iter<maxInter) and (entireSet or (alphaPairsChanged>0))  :
            alphaPairsChanged = 0
            if entireSet :
                for i in range(os.m) :
                    alphaPairsChanged += innerL(i,os)
                iter +=1
            else :
                nonBoundIds = nonzero((os.alphas.A>0)*(os.alphas.A<C))[0]
                for i in nonBoundIds :
                    alphaPairsChanged += innerL(i, os)
                iter +=1
            if entireSet :
                entireSet = False
            elif alphaPairsChanged==0 :
                entireSet = True
        return os.b, os.alphas
    
    dataMat, labelMat = loadDataSet('C:\\2017\\提高\\机器学习\\训练样本\\svm\\testSet.txt')
    b,alphas = smop(dataMat,labelMat,0.6,0.001,40)
    #print(b)
    #print(alphas)
    w = calcWs(alphas, dataMat, labelMat)
    print(w)
    validate(mat(dataMat),labelMat,w,b)

其中最重要的改造是selectJ函数，对j的选择不再是随机。

实测后，改造后的性能提高了15倍左右。

[20180107223139326]: /images/20220602/5d227537ee404af6966793bd6752e83b.png
[20180107223208930]: /images/20220602/005465c8fa8843a49abff55df1e4c38c.png
[20180107223352453]: /images/20220602/9d3ceb4be5484d4e8b8a9e4c5f10b62a.png
[20180107223427728]: /images/20220602/5f67247af6ac4b04a117ab6122c6b909.png
[20180107223459441]: /images/20220602/de17ddbdc80f4b08837fffc37487f79d.png
[20180107223549896]: /images/20220602/a5a6f584104a4a38982b0c0e1a687239.png
[20180107223632192]: /images/20220602/e37c481d9dbf4b77ada0a5ee19df0e64.png
[20180107223658499]: /images/20220602/1d074612d782400ca5d58d6df15aa005.png
[20180107223727940]: /images/20220602/b3e84439703449ffb8bf9acc6c6f5d3a.png
[20180107223803954]: /images/20220602/9037f6478492422a9b3cf0466dd5bc0a.png
[https_www.jianshu.com_p_55458caf0814]: https://www.jianshu.com/p/55458caf0814
[20180107223846840]: /images/20220602/41b44564575a419d87a2a7dc3c35b76c.png
[20180107223916753]: /images/20220602/9f9bf0d77da24256965b4f0bcc3c08a7.png
[20180107223946440]: /images/20220602/d805400c4b28414e8d91530218d7dc35.png
[20180107224015345]: /images/20220602/71a57f6c523a4396b3a6a5d1d8731f13.png
[20180107224059081]: /images/20220602/095eb81d5a9b4221bf930f070d41c79c.png
[20180107224131952]: /images/20220602/0de22c2759e2423689905cb92eb6701a.png
[20180107224202017]: /images/20220602/86470b63c9b94f15a0b20a582a20af92.png
[20180107224227491]: /images/20220602/e79bd39795b94d2c87ac651b3d69ec99.png
[20180107224259751]: /images/20220602/3a9c202e6f2a4a05b206e84cf316c375.png
[20180108132852226]: /images/20220602/5673ca26d7aa48e7bb8b57df05f831f4.png