zl程序教程

您现在的位置是:首页 >  后端

当前栏目

Platt SMO 支持向量机算法(Python实现)

Python算法 实现 支持 向量
2023-09-11 14:22:29 时间

1. SMO SVM算法简述

  1.1 概述

SVM(Support Vector Machine,SVM)算法既是支持适量机算法。算法的原始思想很简单,既是找到一个决策面使两类(本文以两类为例进行说明)分开,且这个决策面和两类之间的间隙尽可能的大。这样带来的好处就是泛化错误率低,能够很好地对两类的问题进行分类。
因而,SVM算法的优点就是泛化错误率低,计算的开销不大,结果容易解释。缺点就是算法对参数的选择和核函数的选择敏感,原始的分类器适合二分类。

  1.2 原理

我们定义决策面为:
点到面的距离:

在SVM算法中距离决策面近的点的g(x)的值为+1(对应类别X1)或是-1(对应类别X2),因而需要满足的条件为:

公式中的X1和X2对应为两个类别,则公式(3)可以被写为:

由于yi有正负对应上式3,则条件全为大于等于1的
因而对于上述的二次优化问题可以有KKT条件,用拉格朗日方程的形式写出:

则将公式(5)中的第一个、第二个式子带入到最后一个式子中,得到:

则对(5)最后一个式子最大化,将(6)带入到(5)中的最后一个式子代数运算之后得到的等价的优化为:

则通过解上诉的方程组,得到决策面的系数,就可以作为以后的分类使用了。但是上述的方法具有缺陷,它是对应下面分布的情况的,在这样的情况中,能够直接找到一个决策面。但是对于接下来的一种分布情况就不适用了(图片来源于度娘)

  1.2 情况2

先来看一下这里分析的分类情况(图片来源于网络):

在这样的分布情况下就出现了3中可能的情况,(1)位于分离段以外并且正确分类的数据;(2)落在分离段内部,但是正确分类了的数据;(3)被错误分类的数据。因而原来的原始约束条件就写为了:

因而,对于情况(1)这里定义的偏移量为0了;情况(2)偏移量就大于0且小于1;情况(3)偏移量大于0。这里引入了偏移量(松弛变量)增加了变量的个数,但是SVM算法的目的是使得分类间隔尽可能大,同时保持偏移量大于0的数据尽可能少。 则原来的最小化函数J就变成了:

其中I(x)是0-1的阶跃函数,但是区间为大于等于0。引入了偏移量之后的推导也是跟之前情况的推导是类似的,则最后得到的优化函数为:


2. Platt的SMO算法

  2.1 SMO方法

SMO方法是将原来SVM的二次规划问题转化为固定大小的二次规划子问题。因为SVM算法需要遵循公式(10)中限制条件1,所以SMO方法选择两个λ进行优化。在SMO方法中使用启发式的方法选择两个λ,这样选取会极大加速算法的收敛速度,做过对比证明要比随机的撞天婚方法快很多很多。这是论文中SMO方法更新λ的部分(偷下懒截出来仅供交流),在该论文的2.2节介绍了循环的条件,读者可以下载论文下来进行阅读。

  2.2 Python实现代码

产生测试数据:
#-*- coding: UTF-8 -*-
import numpy as np
import matplotlib.pyplot as plt

def GenerateData(nums, Distribution_Type = "Guassion"):
    DataCreated = []
    DataLabel = []
    for i in range(nums/2):
        if "Liner" == Distribution_Type:
            temp = 10*np.random.random_sample((2, 1))
        if "Guassion" == Distribution_Type:
            temp = 10*np.random.normal(0, 0.5, (2, 1))
        temp = [temp[0][0], temp[1][0]]
        DataLabel.append(1)
        DataCreated.append(temp)

    for i in range(nums/2):
        if "Liner" == Distribution_Type:
            temp = 10*np.random.random_sample((2, 1))
        if "Guassion" == Distribution_Type:
            temp = 10*np.random.normal(1, 0.5, (2, 1))
        temp = [temp[0][0], temp[1][0]]
        DataLabel.append(-1)
        DataCreated.append(temp)

    DataCreated = np.array(DataCreated)
    DataLabel = np.array(DataLabel)
    #print DataCreated
    #print DataLabel
    ShowData(DataCreated, DataLabel)
    return DataCreated, DataLabel

def ShowData(Data, Datalabel, IsShow = True):
    Data_rows = Data.shape[0]
    Datalabel_rows = Datalabel.shape[0]
    if Data_rows != Datalabel_rows:
        print "ShowData: Data_rows!=Datalabel_rows"
        return
    if IsShow:
        plt.figure("the input data")
        for i in range(Data_rows):
            if 1 == Datalabel[i]:
                plt.plot(Data[i][0], Data[i][1], 'r^')
            if -1 == Datalabel[i]:
                plt.plot(Data[i][0], Data[i][1], 'gs')
        #plt.ion()
        plt.xlabel("X label")
        plt.ylabel("Y label")
        plt.grid(True)
        plt.title("input data distribution")
        plt.show()
    else:
        print "Show Data has been denied"

def SaveData2Txt(Data, Datalabel, filepath = "Data.txt"):
    if "Data.txt" == filepath:
        print "there is no input file-path name, will generate in locally\n"
    try:
        myfile = open(filepath, 'w')
        rows = Data.shape[0]
        for i in range(rows):
            myfile.writelines(Data[i][0].__str__() + "\t" + Data[i][1].__str__() + "\t" +
                              Datalabel[i].__str__() + "\n")
        myfile.close()
    except IOError, e:
        print "write data to %s failed!" % filepath

def ReadData2Mat(filepath = "Data.txt"):
    if "Data.txt" == filepath:
        print "will open %s as Data input\n" % (filepath)
    Data = []
    Label = []
    try:
        myfile = open(filepath, 'r')
        count = 0
        for line in myfile.readlines():
            lineArr = line.strip().split('\t')
            Data.append([float(lineArr[0]), float(lineArr[1])])
            Label.append([float(lineArr[2])])
            count += 1
        myfile.close()
    except IOError, e:
        print "open and read %s error\n" % (filepath)

    Data = np.mat(Data)
    Label = np.mat(Label)
    return Data, Label
SMO方法优化的SVM算法
#-*- coding: UTF-8 -*-
import numpy as np
import matplotlib.pyplot as plt

'''
import functions from other file
'''
import SVM_DataCreate as data_create

'''==============================================================================================
完整的SMO SVM算法
'''

class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = dataMatIn.shape[0]
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m, 2)))

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

'''
函數功能:將輸入的元素限定在一個範圍內
'''
def clipAlpha(input, Low, high):
    if input>high:
        input = high
    if input<Low:
        input = Low

    return input

'''
函數功能:在輸入的參數i和m之間隨機選擇一個不同於i的數字,也就是在選定了i之後隨機選取一個與之配對的alpha的取值的下標
'''
def selectJrand(i, m):
    j = i
    while j==i:
        j = int(np.random.uniform(0, m))

    return j

'''
函數功能:選擇一個SMO算法中與外層配對的alpha值的下標
'''
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:  #啓發式選取配對的j,計算誤差
        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]

'''
SMO算法中的優化部分
'''
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]>oS.C)):
        j, Ej = selectJ(i, oS, Ei)
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        '''公式(7)'''
        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 H == L:
            print "H==L program continued"
            return 0

        '''公式(8)(9)'''
        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 0 <= eta:
            print "eta>=0 program continued"
            return
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], L, H)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print "j not moving enough %s" % ("program continued")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        updateEk(oS, i) #更新誤差緩存'

        '''設置常數項 b '''
        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

'''
完整版SMO算法
dataMatIn: 訓練數據
classLabels: 數據標籤
C: 常量
toler: 容錯度
maxIter: 最大迭代次數
kTup=('lin', 0): 核函數類型
'''
def SMOP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(np.mat(dataMatIn), np.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 = np.nonzero((oS.alphas.A>0) * (oS.alphas.A<C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print "fullSet, iter: %d    i:%d, pairs changed: %d" % (iter, i, alphaPairsChanged)
            iter += 1
        if entireSet:
            entireSet = False
        elif (0==alphaPairsChanged):
            entireSet = True
        print "iteration number: %d" % (iter)
    return oS.b, oS.alphas


'''
函數功能:由計算出來的alphas獲得進行分類的權重向量
alphas: 計算出來的alpha向量
dataArr: 訓練數據
classLabels: 數據標籤
'''
def calcWs(alphas, dataArr, classLabels):
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels)
    rows, cols = np.shape(dataArr)
    w = np.mat(np.zeros((cols, 1)))
    for i in range(rows):
        w += np.multiply(alphas[i, :]*labelMat[:, i], X[i,:].T)
    return w

'''
main function
'''
Data, DataLabel = data_create.GenerateData(10)     #generate data
print Data
b, alphas = SMOP(Data, DataLabel, 0.6, 0.001, 50)

SuportPoint = []
for i in range(np.shape(alphas)[0]):
    if 0 < alphas[i]:
        SuportPoint.append([i, alphas[i]])
SuportPoint = np.mat(np.array(SuportPoint))
print SuportPoint
weight_vector = calcWs(alphas, np.mat(Data), np.mat(DataLabel))
print weight_vector

#fig, ax =plt.subplots(1,2,figsize=(10,5))
plt.figure("decsion boundry")
plt.title("decsion boundry")
p1 = plt.subplot(121)
p2 = plt.subplot(122)
for i in range(Data.shape[0]):
    if 1 == DataLabel[i]:
        p1.plot(Data[i][0], Data[i][1], 'r^')
    if -1 == DataLabel[i]:
        p1.plot(Data[i][0], Data[i][1], 'gs')
p1.set_xlabel("X label")
p1.set_ylabel("Y label")
x = np.arange(-10, 10, 1)
l = float( weight_vector[0, :])
k = np.multiply(x, l)
y = np.multiply(weight_vector[0, :], x) + weight_vector[1, :] + b
y = np.array(y)
p1.plot(x, y[0, :], 'b')
p1.grid(True)

x = np.arange(0, 10, 0.5)
l = float( weight_vector[0, :])
k = np.multiply(x, l)
y = np.multiply(weight_vector[0, :], x) + weight_vector[1, :] + b
y = np.array(y)
p2.plot(x, y[0, :], 'b')
p2.set_xlabel("X label")
p2.set_ylabel("Y label")
p2.grid(True)
plt.title("input data distribution")
'''
SaveImg = raw_input("保存圖片(0), 不保存圖片(1)")
if "0" == SaveImg:
    try:
        plt.savefig('result.png', format='png')
    except Exception, e:
        print "save Image Error: %s" % (e)
elif "1" == SaveImg:
    print "result image has not been saved"
else:
    print "Input Error, please input '1' or '0' "
'''
plt.show()

print b
print alphas


运行结果


3. 核函数方法

当遇到下面这种情况(线性不可分割)的时候就需要引入核函数的概念(图片仅作为说明数据的分布情况,来源于网络)

核函数说白了就是一个空间的映射,从一个特征空间映射到另外一个特征空间,一般情况下是从低维度到高维度的。将特征进行映射操作之后,原本线性不可分的数据,在高维的空间中变成了可以线性可分的了。在一般情况下使用的比较多的是径向基核函数,此外还有一些核函数,如傅里叶核等,读者可查询其它相关资料

那么引入核函数之后,对于上面的代码需要更改部分
添加函数
'''==============================================================================================
完整的SMO SVM算法
'''
'''函數說明:核函數方法
DataIn: 訓練數據集合輸入
Sample: 輸入與的樣本數據
kTup: 核函數的類型
'''
def KernelTrans(DataIn, Sample, kTup):
    DataIn = np.mat(DataIn)
    rows, cols = np.shape(DataIn)
    k = np.zeros((rows, 1))
    if "lin" == kTup[0]:
        k = DataIn * Sample.T
    elif "rbf" == kTup[0]:
        for i in range(rows):
            deltaRow = DataIn[i, :] - Sample
            k[i] = deltaRow*deltaRow.T
        k = np.exp(k / (-1*kTup[1]**2))
    else:
        print "not legal input kernel"
    return k

修改:
'''SMOP算法算法結構體定義'''
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup = ["rbf", 1]):
<span style="white-space:pre">	</span>...
        self.K = np.mat(np.zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = KernelTrans(self.X, self.X[i, :], kTup)

        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
</pre><pre name="code" class="python">        '''設置常數項 b '''
        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]

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

4. 参考资料

1. 《机器学习实战》
2. 《模式识别》 Sergios Theodoridis