SVM SMO算法代码详细剖析

SVM SMO算法代码详细剖析算法实现一:本文要结合SVM理论部分来看即笔者另一篇https://blog.csdn.net/weixin_42001089/article/details/83276714二:有了理论部分下面就是直接代码啦,本文用四部分进行介绍:最简版的SMO,改进版plattSMO,核函数,sklearn库的SVM,四部分以%%%%%%%分开,采取的顺序是先给代码及结果,然后分析三:这里代码大…

大家好,又见面了,我是你们的朋友全栈君。

算法实现

一:本文要结合SVM理论部分来看即笔者另一篇svm原理详细推导_爱吃火锅的博客-CSDN博客_svm推导

二:有了理论部分下面就是直接代码啦,本文用四部分进行介绍:最简版的SMO,改进版platt SMO,核函数,sklearn库的SVM,四部分以%%%%%%%分开,采取的顺序是先给代码及结果,然后分析

三:这里代码大部分来自于Peter Harrington编写的Machine Learning in Action

其网络资源Manning | Machine Learning in Action

四:代码中需要注意的一点就是采用启发式来寻找需要优化的\alpha

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

首先看一下最简版的smo优化算法

这里有两个py文件,一个是用来构造SVM的,一个是用来测试的:

MySVM:

# -*- coding: utf-8 -*-
import random
import  numpy  as np
import matplotlib.pyplot as plt
#辅助函数一
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:
        aj = H
    if L > aj:
        aj = L
    return aj

#最简版本SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = np.mat(dataMatIn); 
    labelMat = np.mat(classLabels).transpose()
    b = 0; 
    m,n = np.shape(dataMatrix)
    alphas = np.mat(np.zeros((m,1)))
    iter_num = 0
    while (iter_num < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            #注意一
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[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)):
                j = selectJrand(i,m)
                fXj = float(np.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]):
                    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 - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: 
                    print("eta>=0"); 
                    continue
                
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
              
                alphas[j] = clipAlpha(alphas[j],H,L)
               
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化小,不需要更新"); continue
                #注意三
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])

                #注意四
                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
              
                print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
       
        if (alphaPairsChanged == 0): 
            iter_num += 1
        else: iter_num = 0
        print("迭代次数: %d" % iter_num)
    #注意五
    return b,alphas
def calcWs(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
    w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()

def showClassifer(dataMat,labelMat,alphas, w, b):
    data_plus = []                                  
    data_minus = []
    #注意六
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)              
    data_minus_np = np.array(data_minus)            
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)   
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) 
    #注意七
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #注意八
    for i, alpha in enumerate(alphas):
        if 0.6>abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
        if 0.6==abs(alpha) :
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='yellow')
    plt.show()

接着是测试函数(MyTest):

# -*- coding: utf-8 -*-#
import MySVM as svm
x=[[1,8],[3,20],[1,15],[3,35],[5,35],[4,40],[7,80],[6,49],[1.5,25],[3.5,45],[4.5,50],[6.5,15],[5.5,20],[5.8,74],[2.5,5]]
y=[1,1,-1,-1,1,-1,-1,1,-1,-1,-1,1,1,-1,1]
b,alphas = svm.smoSimple(x,y,0.6,0.001,40)
w = svm.calcWs(x,y,alphas)
svm.showClassifer(x,y,alphas, w, b)

运行结果:

SVM SMO算法代码详细剖析

,,,,,,,,,,,,,

SVM SMO算法代码详细剖析

SVM SMO算法代码详细剖析

接下来一步步分析MySVM中的代码

首次看两个简单的辅助函数,第一个函数的作用就是用来选择\alpha对的(即寻找i,j这一对)

第二个函数就是为了将\alpha规划到[0,C]范围内,对应到理论推导部分的:

w=\sum_{i=1}^{m}\alpha _{i}y_{i}x_{i}\, \, \, \, (9)

可能这里和上篇最后给出w的更新形式上看上去有点不对应,其实是一样的,推导部分即最后一张图是:

w = w+y_{1}(\alpha_{1} ^{*}-\alpha_{1} )x_{1}+y_{2}(\alpha_{2} ^{*}-\alpha_{2} )x_{2}

而这里的其实在程序就是对应的就是如下三步:

 fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b

 alphas[j] -= labelMat[j]*(Ei - Ej)/eta
 alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
 

紧接着的if这里就是启发式选择,即寻找那些误差过大(正间隔和否间隔)且在(0,C)范围内的\alpha进行优化,选择误差大的进行优化我们很容易理解,那为什么要选择(0,C)范围内,而不选择边界值呢(值等于0或C),那是因为它们已经在边界啦,因此不再能够减少或者增大具体细节请看推导部分,该部分包括L和H为什么要这样赋值,以及为什么L==H的时候要返回都有讲到。

注意二的部分就是类似我们在推导部分的\alpha _{2}更新,只不过这里有一步如果变化太小我们就不更新了,直接跳过,只不过原公式中的\theta =-eta,所以推导中原本\alpha _{2}是加上后面的,而这里是减即:本质是一样的啦,当然\alpha _{1}更新方向要相反,所以代码中对应的是+

 alphas[j] -= labelMat[j]*(Ei - Ej)/eta

至于为什么eta >= 0为什么要跳过该次循环,请看推导部分,只不过因为\theta =-eta所以原来是过滤掉<=0,这里是>=0

注意三部分就是类似\alpha _{1}的更新即大小和\alpha _{2}相同,方向相反

注意四的部分应该很直观啦,看我们推导的b的更新结论一目了然

注意五返回的b是一个实数,alphas是一个[m,1]矩阵

最后说一下smo函数的alphaPairsChanged和iter_num以及maxIter的参数意义,maxIter是最外部的大循环,是人为设定的最大循环次数,循环为最大次数后就强行结束返回,在每一个大循环下都有一个for循环,用以遍历一遍所有的\alpha,遍历完这一遍所有的\alpha后,alphaPairsChange用以记录看有多少对被优化啦,如果alphaPairsChange不为0,即这一遍走下来后,我们进行了优化,也就代表目前\alpha还不够好,所以我们将iter_num设为0,继续优化,当alphaPairsChange为0时,说明我们这一遍走下来,说明\alpha都很好啦,没有优化的必要啦,我们将iter_num加一,接着下一遍再去整体看看\alpha,如果还是alphaPairsChange为0,恩恩,不错,不错,将iter_num再加一,如果iter_num到了maxIter即连续进行了maxIter遍整体(for)观察都没发现需要优化的\alpha,说明足够好了,返回吧!!!!!!一旦中间出现意外,发现有需要优化的,就至少说明有不完善的地方,那么我们立马让iter_num为0,即一定要达到连续遍历maxIter次都没发现不足,我们才放心,才返回,发现瑕疵立马iter_num=0从头开始,怎么样?就是这么严格,这也是上面运行结果开始的时候为什么迭代次数一直都是0,后面趋于收敛,迭代次数连续增加,直到maxIter结束返回

正是因为如此,可以想象得的到带来的结果就是 时间复杂度太高,所以有了后来改进版本的Platt SMO,后面介绍

接下来的calcWs函数作用是:根据训练出来的\alpha生成w:

对应的公式就是:

w=\sum_{i=1}^{m}\alpha _{i}y_{i}x_{i}\, \, \, \, (9)

目前我们已经训练除了SVM模型,即得到了我们想要的w和b,对应的步骤就在上面所说的黄色部分

接下来可视化看一下结果showClassifer:

注意六的部分就是我们把原始点画出来,不同的颜色代表不同的分离(橙色的点对应的标签是-1,蓝色的点对应的标签是1)

注意七的部分就是画出训练出来的超平面

y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2

这个很好理解啦,超平面是:

0=x_{1}w_{1}+x_{2}w_{2}+b

所以:

x_{2}=\frac{-x_{1}w_{1}-b}{w_{2}}

程序中为了让超平面尽可能的横穿整个数据点,所以选取了所有点中x坐标最大和最小的点即x1和x2:

然后利用上面的公式,计算出了对应的纵坐标

注意八的部分就是画出向量机点,和那些我们“忽略”的点,依据是推导的:

\begin{matrix} \alpha _{i}=0\Leftrightarrow y_{i}y(x)\geq 1\\ 0< \alpha _{i}< C\Leftrightarrow y_{i}y(x)= 1 \\ \alpha _{i}=C\Leftrightarrow y_{i}y(x)\leq 1 \end{matrix}

即在点在两条间隔线外则,对应前面的系数\alpha为0,在两条间隔线里面的对应的系数为C,在两条间隔线上的点对应的系数在0和C之间。至于为什么请看上篇的推导细节

带有红色圆圈的是支持向量机点即间隔线上的点,带有黄色的点是间隔线内的点

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Platt SMO

其是SMO算法的一个改进版,速度更快。

其主要变化的地方有两个

一:在使用启发式方法选择了一个\alpha后,我们会去选择另外一个与之对应是吧,即

 j = selectJrand(i,m)

但是改进的的SMO算法中,这里也使用启发式来选择,即选择与Ei误差最大的Ej即选择最大步长,简单来说就是找最需要优化的j,而不是像最简版本那样,毫无目的的随机去选择,所以对应到推导公式里面就是\alpha _{1}\alpha _{2}都采用启发式来寻找

二:改进后的算法是采用在“非边界值”和“边界值”范围内交替遍历优化的

下面来看一下具体代码:

smoP:

# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
import random
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 
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  
        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)))
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:
        aj=H
    if L>aj:
        aj=L
    return aj
        
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):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1,Ei]
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]  
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE): 
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej
def updateEk(oS, k): 
    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)
        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
        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)
        if (abs(oS.alphas[j] - alphaJold) < oS.tol): 
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        updateEk(oS, i) 
        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]<oS.C):
            oS.b = b1
        elif (0 < oS.alphas[j]<oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0

def calcWs(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = array(alphas), array(dataMat), array(labelMat)
    w = dot((tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): 
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    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
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas

def showClassifer(dataMat,labelMat,alphas, w, b):
    data_plus = []                                  
    data_minus = []
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = array(data_plus)              
    data_minus_np = array(data_minus)            
    plt.scatter(transpose(data_plus_np)[0], transpose(data_plus_np)[1], s=30, alpha=0.7)   
    plt.scatter(transpose(data_minus_np)[0], transpose(data_minus_np)[1], s=30, alpha=0.7) 
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    for i, alpha in enumerate(alphas):
        if 0.6>abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
        if 50==abs(alpha) :
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='yellow')
    plt.show()


接着还是测试函数(MyTest):

# -*- coding: utf-8 -*-#
import smoP as svm
x=[[1,8],[3,20],[1,15],[3,35],[5,35],[4,40],[7,80],[6,49],[1.5,25],[3.5,45],[4.5,50],[6.5,15],[5.5,20],[5.8,74],[2.5,5]]
y=[1,1,-1,-1,1,-1,-1,1,-1,-1,-1,1,1,-1,1]
b,alphas = svm.smoP(x,y,50,0.001,40)
w = svm.calcWs(x,y,alphas)
svm.showClassifer(x,y,alphas, w, b)

SVM SMO算法代码详细剖析

,,,,,,,,,,,

SVM SMO算法代码详细剖析

SVM SMO算法代码详细剖析

这里首先optStruct函数定义了一个类作为数据结构来存储一些信息,这里面的alphas就是我们的\alpha,eCache第一列就是一个是否有效的标志位,第二列存储着误差值E,总之这个结构体的定义就是为了作为一个整体,方便调用,管理。

calcEk和最简版本没什么差别,只不过我们已经定义了结构体,所以直接可以调用结构体便可得到一些信息,所以下面所有代码都是这样,比如C我们可以直接用oS.C等等

selectJ和最简版本不一样啦,这里也就是我们说的用启发式来寻找j,这里的:

if (len(validEcacheList)) > 1:

主要是防止第一次循环的时候,如果是第一次那么就随机选择,之后都使用启发式来选择

 updateEk就是用来在计算完i和j的Ei和Ej后更新数据结构中的的eCache

innerL和最简版本的smoSimple内循环(就是for循环下面的代码:用来优化\alpha和b的核心代码)一模一样,只不过这里要把一些东西,改为数据结构中定义的,而且这里的selectJ已经采用了启发式寻找

接下来就是我们的smoP,也是Platt SMO利用主循环封装整个算法的过程,其和最简版本一样,也是两个循环:

外训练也是使用了一个maxIter,同时使用了iter来记录遍历次数(对应最简版本的iter_num),但是两者含义却不一样,这里的iter就是单纯的代表一次循环,而不管循环内部具体做了什么,它没有被清0这个过程,随着程序运行一直是个累加的过程,上面运行结果也可以看到iter是一直递增的,这也是Platt SMO之所以能够加快算法的一个重要原因,而最简版本的iter_num要肩负着连续这一条件,同时这里的外循环相对于最简版本的的外循环多了一个退出条件即:遍历整个集合都没发现需要改变的\alpha(说明都优化好啦,退出吧)

再来看一下内循环,这里对应着两种情况,一种是在全集上面遍历([0,C]),另一种是非边界上面((0,C)),通过

if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True

使两种情况交替遍历

其他部分包括W的获得,可视化什么的就和最简版本一样啦,不再重复介绍啦

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

核函数

核函数的作用细节请看推导部分,核函数种类很多,这里看一下最常用的径向基高斯(RBF)核函数 

下面来简单说一下部分代码(这里只说不同的地方,相同的地方不再重述)

# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
def loadDataSet(filename): 
    dataMat=[]
    labelMat=[]
    fr=open(filename)
    for line in fr.readlines():
        lineArr=line.strip().split(',')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat 

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:
        aj=H
    if L>aj:
        aj=L
    return aj

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


class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  
        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))) 
        self.K = mat(zeros((self.m,self.m))) 
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)


def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1,Ei]
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]  
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE): 
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


def updateEk(oS, k): 
    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) 
        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.K[i,j] - oS.K[i,i] - oS.K[j,j] 
        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)
        if (abs(oS.alphas[j] - alphaJold) < oS.tol): 
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        updateEk(oS, i) 
       
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
        if (0 < oS.alphas[i]<oS.C):
            oS.b = b1
        elif (0 < oS.alphas[j]<oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0


def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    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
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas


def testRbf(data_train,data_test):
    dataArr,labelArr = loadDataSet(data_train) 
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 0.2)) 
    datMat=mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas)[0]  
    sVs=datMat[svInd] 
    labelSV = labelMat[svInd] 
    print("there are %d Support Vectors" % shape(sVs)[0]) 
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', 1.3)) 
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b  
        if sign(predict)!=sign(labelArr[i]): 
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m)) 
    dataArr_test,labelArr_test = loadDataSet(data_test) 
    errorCount_test = 0
    datMat_test=mat(dataArr_test)
    labelMat = mat(labelArr_test).transpose()
    m,n = shape(datMat_test)
    for i in range(m): 
        kernelEval = kernelTrans(sVs,datMat_test[i,:],('rbf', 0.1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr_test[i]):
            errorCount_test += 1
    print("the test error rate is: %f" % (float(errorCount_test)/m))
    
    return dataArr,labelArr,alphas

def showClassifer(dataMat,labelMat,alphas):
    data_plus = []                                  
    data_minus = []
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = array(data_plus)              
    data_minus_np = array(data_minus)            
    plt.scatter(transpose(data_plus_np)[0], transpose(data_plus_np)[1], s=30, alpha=0.7)   
    plt.scatter(transpose(data_minus_np)[0], transpose(data_minus_np)[1], s=30, alpha=0.7) 
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()



    

MyTest:

# -*- coding: utf-8 -*-
import smoPrbf as svm
traindata='C:\\Users\\asus-\\Desktop\\train_data.csv'
testdata='C:\\Users\\asus-\\Desktop\\test_data.csv'
TraindataArr,TrainlabelArr,alphas = svm.testRbf(traindata,testdata)
svm.showClassifer(TraindataArr,TrainlabelArr,alphas)

\delta =0.2时:

SVM SMO算法代码详细剖析

SVM SMO算法代码详细剖析

 当\delta =1.2时:

SVM SMO算法代码详细剖析

SVM SMO算法代码详细剖析

kernelTrans函数的作用就是核函数的计算部分,对应到推导公式是:

K(x_{i},x)=exp(-\frac{\left \| x_{i}-x \right \|^{2}}{2\delta ^{2}})

这里的kTup就是指定使用什么核函数,kTup[0]参数是核函数类型,kTup[1]是核函数需要的超参数,注意这里只支持线性和径向基高斯(RBF)核函数 两种.

optStruct函数增加了一个字段即K,其是一个m*m的矩阵。注意它的含义:

我们的核函数是K(x_{i},x)即拿来一个点x,要和所有样本x_{i}做运算,这里的行代表的意义就是所有样本,列代表的是x所以这里是:

elf.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

innerL变化的部分是:

eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]

对比之前的:

 eta = 2.0 * oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T-oS.X[j,:]*oS.X[j,:].T

就可以更好的理解为什么在optStruct结构体中K字段要这么设计

同理变化的地方还有:

 b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
 b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]

calcEk变化的地方:

fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)

简单来说就是原先有x_{i}x的地方都要换成核函数的内积形式即

注意这里在通过\alpha构建权重w时是只用到是支持向量机那些点即

程序也会将向量机的个数打印出来。

最后讨论一下rbf的超参数意义即\delta值对应在在代码的:

b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 0.2)) 

b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 1.2)) 

通过上面的实验我们可以大体看出随着\delta的增加,支持向量机的个数在减少,由原来的43个减少到25(红色圆圈的点就是支持向量机),\delta的取值存在一个最优解,当\delta太大,支持向量机太少,也就是说我们利用了很少的点去决策,显然结果不好,正如上面体现的那样,测试集的错误在上升,当\delta太小,支持向量机太多,也就是我们基本利用了所有样本点,其实这个时候已经退化到类似KNN啦,因为KNN就是利用了到所有样本点的距离来决策的,可能会有这样的疑问?不对呀?使用KNN时不是会指定利用多少个点吗?不是利用所有点呀?哈哈哈,仔细想想,它的过程是先计算到全部样本的距离,然后再从中选择K个最近距离的点来进行比较的,所以它每次要用到的是全部样本点,而SVM是一旦训练出\alpha后,在之后的决策中就只使用

关于更多sklearn的SVM更多调用请看:

python_sklearn机器学习算法系列之SVM支持向量机算法_爱吃火锅的博客-CSDN博客

当然还有其他机器学习的调用

# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
from sklearn import metrics
from sklearn import svm

def loadDataSet(filename): 
    dataMat=[]
    labelMat=[]
    fr=open(filename)
    for line in fr.readlines():
        lineArr=line.strip().split(',')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat


traindata='C:\\Users\\asus-\\Desktop\\train_data.csv'
testdata='C:\\Users\\asus-\\Desktop\\test_data.csv'


x_train,y_train = loadDataSet(traindata)
x_test,y_test = loadDataSet(testdata)


clf1 = svm.SVC(C=0.8, kernel='rbf', gamma=10, decision_function_shape='ovr')
clf2 = svm.SVC(C=0.8, kernel='rbf', gamma=20, decision_function_shape='ovr')

clf1.fit(x_train, y_train)
clf2.fit(x_train, y_train)



y_predict1=clf1.predict(x_test)
y_predict2=clf2.predict(x_test)
print('--------------------- gamma=10--------------------------------')
print(metrics.classification_report(y_test,y_predict1))
print('--------------------- gamma=20--------------------------------')
print(metrics.classification_report(y_test,y_predict2))

SVM SMO算法代码详细剖析

结尾给一下我们用的数据集,方便大家实验:

Train_data:

-0.214824 0.662756 -1
-0.061569 -0.091875 1
0.406933 0.648055 -1
0.22365 0.130142 1
0.231317 0.766906 -1
-0.7488 -0.531637 -1
-0.557789 0.375797 -1
0.207123 -0.019463 1
0.286462 0.71947 -1
0.1953 -0.179039 1
-0.152696 -0.15303 1
0.384471 0.653336 -1
-0.11728 -0.153217 1
-0.238076 0.000583 1
-0.413576 0.145681 1
0.490767 -0.680029 -1
0.199894 -0.199381 1
-0.356048 0.53796 -1
-0.392868 -0.125261 1
0.353588 -0.070617 1
0.020984 0.92572 -1
-0.475167 -0.346247 -1
0.074952 0.042783 1
0.394164 -0.058217 1
0.663418 0.436525 -1
0.402158 0.577744 -1
-0.449349 -0.038074 1
0.61908 -0.088188 -1
0.268066 -0.071621 1
-0.015165 0.359326 1
0.539368 -0.374972 -1
-0.319153 0.629673 -1
0.694424 0.64118 -1
0.079522 0.193198 1
0.253289 -0.285861 1
-0.035558 -0.010086 1
-0.403483 0.474466 -1
-0.034312 0.995685 -1
-0.590657 0.438051 -1
-0.098871 -0.023953 1
-0.250001 0.141621 1
-0.012998 0.525985 -1
0.153738 0.491531 -1
0.388215 -0.656567 -1
0.049008 0.013499 1
0.068286 0.392741 1
0.7478 -0.06663 -1
0.004621 -0.042932 1
-0.7016 0.190983 -1
0.055413 -0.02438 1
0.035398 -0.333682 1
0.211795 0.024689 1
-0.045677 0.172907 1
0.595222 0.20957 -1
0.229465 0.250409 1
-0.089293 0.068198 1
0.3843 -0.17657 1
0.834912 -0.110321 -1
-0.307768 0.503038 -1
-0.777063 -0.348066 -1
0.01739 0.152441 1
-0.293382 -0.139778 1
-0.203272 0.286855 1
0.957812 -0.152444 -1
0.004609 -0.070617 1
-0.755431 0.096711 -1
-0.526487 0.547282 -1
-0.246873 0.833713 -1
0.185639 -0.066162 1
0.851934 0.456603 -1
-0.827912 0.117122 -1
0.233512 -0.106274 1
0.583671 -0.709033 -1
-0.487023 0.62514 -1
-0.448939 0.176725 1
0.155907 -0.166371 1
0.334204 0.381237 -1
0.081536 -0.106212 1
0.227222 0.527437 -1
0.75929 0.33072 -1
0.204177 -0.023516 1
0.577939 0.403784 -1
-0.568534 0.442948 -1
-0.01152 0.021165 1
0.87572 0.422476 -1
0.297885 -0.632874 -1
-0.015821 0.031226 1
0.541359 -0.205969 -1
-0.689946 -0.508674 -1
-0.343049 0.841653 -1
0.523902 -0.436156 -1
0.249281 -0.71184 -1
0.193449 0.574598 -1
-0.257542 -0.753885 -1
-0.021605 0.15808 1
0.601559 -0.727041 -1
-0.791603 0.095651 -1
-0.908298 -0.053376 -1
0.12202 0.850966 -1
-0.725568 -0.292022 -1

Test_data:

0.676771 -0.486687 -1
0.008473 0.18607 1
-0.727789 0.594062 -1
0.112367 0.287852 1
0.383633 -0.038068 1
-0.927138 -0.032633 -1
-0.842803 -0.423115 -1
-0.003677 -0.367338 1
0.443211 -0.698469 -1
-0.473835 0.005233 1
0.616741 0.590841 -1
0.557463 -0.373461 -1
-0.498535 -0.223231 -1
-0.246744 0.276413 1
-0.76198 -0.244188 -1
0.641594 -0.479861 -1
-0.65914 0.52983 -1
-0.054873 -0.2389 1
-0.089644 -0.244683 1
-0.431576 -0.481538 -1
-0.099535 0.728679 -1
-0.188428 0.156443 1
0.267051 0.318101 1
0.222114 -0.528887 -1
0.030369 0.113317 1
0.392321 0.026089 1
0.298871 -0.915427 -1
-0.034581 -0.133887 1
0.405956 0.20698 1
0.144902 -0.605762 -1
0.274362 -0.401338 1
0.397998 -0.780144 -1
0.037863 0.155137 1
-0.010363 -0.00417 1
0.506519 0.486619 -1
0.000082 -0.020625 1
0.057761 -0.15514 1
0.027748 -0.553763 -1
-0.413363 -0.74683 -1
0.0815 -0.014264 1
0.047137 -0.491271 1
-0.267459 0.02477 1
-0.148288 -0.532471 -1
-0.225559 -0.201622 1
0.77236 -0.518986 -1
-0.44067 0.688739 -1
0.329064 -0.095349 1
0.97017 -0.010671 -1
-0.689447 -0.318722 -1
-0.465493 -0.227468 -1
-0.04937 0.405711 1
-0.166117 0.274807 1
0.054483 0.012643 1
0.021389 0.076125 1
-0.104404 -0.914042 -1
0.294487 0.440886 -1
0.107915 -0.493703 -1
0.076311 0.43886 1
0.370593 -0.728737 -1
0.40989 0.306851 -1
0.285445 0.474399 -1
-0.870134 -0.161685 -1
-0.654144 -0.675129 -1
0.285278 -0.76731 -1
0.049548 -0.000907 1
0.030014 -0.093265 1
-0.128859 0.278865 1
0.307463 0.085667 1
0.02344 0.298638 1
0.05392 0.235344 1
0.059675 0.533339 -1
0.817125 0.016536 -1
-0.108771 0.477254 1
-0.118106 0.017284 1
0.288339 0.195457 1
0.567309 -0.200203 -1
-0.202446 0.409387 1
-0.330769 -0.240797 1
-0.422377 0.480683 -1
-0.295269 0.326017 1
0.261132 0.046478 1
-0.492244 -0.319998 -1
-0.384419 0.09917 1
0.101882 -0.781145 -1
0.234592 -0.383446 1
-0.020478 -0.901833 -1
0.328449 0.186633 1
-0.150059 -0.409158 1
-0.155876 -0.843413 -1
-0.098134 -0.136786 1
0.110575 -0.197205 1
0.219021 0.054347 1
0.030152 0.251682 1
0.033447 -0.122824 1
-0.686225 -0.020779 -1
-0.911211 -0.262011 -1
0.572557 0.377526 -1
-0.073647 -0.519163 -1
-0.28183 -0.797236 -1
-0.555263 0.126232 -1

看到很多小伙伴私信和关注,为了不迷路,欢迎大家关注笔者的微信公众号,会定期发一些关于NLP的干活总结和实践心得,当然别的方向也会发,一起学习:

SVM SMO算法代码详细剖析
​​​​​​​

 

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/149631.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号