svmMLiA.py据某，為沒有用啟發(fā)式算法峦失，隨機(jī)選擇alphas[i],alphas[j]的SMO算法的實(shí)現(xiàn)征绎。
svmQuicken.py，為啟用了啟發(fā)是算法選擇alphas[i],alphas[j]的SMO算法的實(shí)現(xiàn)提茁。
代碼寫的有點(diǎn)亂，結(jié)果出來之前馁菜，沒心思整理代碼茴扁，結(jié)果出來后，就更沒心思整理代碼了汪疮。
（以下正確率的結(jié)果峭火，都是由訓(xùn)練數(shù)據(jù)獲得超平面之后毁习，再拿訓(xùn)練數(shù)據(jù)去測(cè)試的。沒有專門去整理測(cè)試數(shù)據(jù)）
一卖丸、線性可分-線性核
下圖一是svmMLiA.py實(shí)現(xiàn)的蜓洪，線性可分，隨機(jī)選擇alphas[i]和alphas[j]坯苹，很慢隆檀，但效果很好。紅色點(diǎn)大圈為支持向量點(diǎn)粹湃，正確率為0.92
圖一
下圖二是svmQuicken.py實(shí)現(xiàn)的恐仑，線性可分，啟發(fā)式算法選擇alphas[i]和alphas[j]为鳄，很快裳仆，效果還行，不如圖一孤钦。紅色點(diǎn)大圈為支持向量點(diǎn)歧斟，正確率為0.90
圖二
二、線性不可分-高斯核
下圖都是svmMLiA.py實(shí)現(xiàn)的偏形，隨機(jī)選擇alphas[i]和alphas[j]静袖。圖三、圖四俊扭、圖五分別是高斯核的參數(shù)sigma = 0.1队橙、0.3、0.6得出的結(jié)果萨惑，紅色的大圈為支持向量捐康。顯然，sigma越小庸蔼，得到的支持向量點(diǎn)越多解总，結(jié)果越準(zhǔn)確，如果支持向量太多姐仅，相當(dāng)于每次都利用整個(gè)數(shù)據(jù)集進(jìn)行分類花枫，這時(shí)便成了K近鄰算法了。sigma = 0.1萍嬉、0.3乌昔、0.6的準(zhǔn)確率分別是0.915254，0.830508壤追，0.813559
圖三
圖四
圖五
svmMLiA.py
'''
Created on 2017年7月9日
@author: fujianfei
'''
from os.path import os
import numpy as np
#導(dǎo)入數(shù)據(jù)磕道，數(shù)據(jù)集
def loadDataSet (fileName):
data_path = os.getcwd()+'\\data\\'
labelMat = []
svmData = np.loadtxt(data_path+fileName,delimiter=',')
dataMat=svmData[:,0:2]
#零均值化
# meanVal=np.mean(dataMat,axis=0)
# dataMat=dataMat-meanVal
label=svmData[:,2]
for i in range (np.size(label)):
if label[i] == 0 or label[i] == -1:
labelMat.append(float(-1))
if label[i] == 1:
labelMat.append(float(1))
return dataMat.tolist(),labelMat
#簡(jiǎn)化版SMO算法，不啟用啟發(fā)式選擇alpha行冰，先隨機(jī)選
def selectJrand(i,m):
j=i
while (j==i):
j = int(np.random.uniform(0,m))#在0-m的隨機(jī)選一個(gè)數(shù)
return j
#用于調(diào)整大于H或小于L的值溺蕉，剪輯最優(yōu)解
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
'''
定義核函數(shù)
kernelOption=linear 線性
kernelOption=rbf 高斯核函數(shù)
'''
def calcKernelValue(matrix_x, sample_x, kernelOption):
kernelType = kernelOption[0]
numSamples = matrix_x.shape[0]
kernelValue = np.mat(np.zeros((numSamples, 1)))
if kernelType == 'linear':
kernelValue = matrix_x.dot(sample_x.T)
elif kernelType == 'rbf':
sigma = kernelOption[1]
if sigma == 0:
sigma = 1.0
for i in range(numSamples):
diff = matrix_x[i, :] - sample_x
kernelValue[i] = np.exp(diff.dot(diff.T) / (-2.0 * sigma**2))
else:
raise NameError('Not support kernel type! You can use linear or rbf!')
return kernelValue
'''
簡(jiǎn)化版SMO算法伶丐。
dataMatIn:輸入的數(shù)據(jù)集
classLabels:類別標(biāo)簽
C:松弛變量前的常數(shù)
toler:容錯(cuò)率
maxIter:最大循環(huán)數(shù)
'''
def smoSimple(dataMatIn,classLabels,C,toler,maxIter,kernelOption):
dataMatrix = np.mat(dataMatIn);labelMat = np.mat(classLabels).transpose()
b=0;m,n = np.shape(dataMatrix)
alphas = np.mat(np.zeros((m,1)))
iter = 0
while(iter < maxIter):
alphaPairsChanged = 0 #記錄alpha值是否優(yōu)化，即是否變化
for i in range(m):#遍歷數(shù)據(jù)集疯特，第一層循環(huán)哗魂，遍歷所有的alpha
fXi = float(np.multiply(alphas,labelMat).T.dot(calcKernelValue(dataMatrix,dataMatrix[i,:],kernelOption))) + 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.dot(calcKernelValue(dataMatrix,dataMatrix[j,:],kernelOption))) + 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 * calcKernelValue(dataMatrix[i,:],dataMatrix[j,:],kernelOption) - calcKernelValue(dataMatrix[i,:],dataMatrix[i,:],kernelOption) - calcKernelValue(dataMatrix[j,:],dataMatrix[j,:],kernelOption)
if(eta >= 0):print('eta >= 0');('alpha[j]=%f###############################' % alphas[j]);continue
alphas[j] -= labelMat[j] * (Ei - Ej)/eta
alphas[j] = clipAlpha(alphas[j], H, L)
if (abs(alphas[j]-alphaJold) < 0.0001) : print('j not moving enough');continue
alphas[i] += labelMat[i]*labelMat[j]*(alphaJold - alphas[j])
b1 = b - Ei - labelMat[i]*(alphas[i] - alphaIold)*calcKernelValue(dataMatrix[i,:],dataMatrix[i,:],kernelOption) - labelMat[j]*(alphas[j]-alphaJold)*calcKernelValue(dataMatrix[j,:],dataMatrix[i,:],kernelOption)
b2 = b - Ej - labelMat[i]*(alphas[i] - alphaIold)*calcKernelValue(dataMatrix[i,:],dataMatrix[j,:],kernelOption) - labelMat[j]*(alphas[j]-alphaJold)*calcKernelValue(dataMatrix[j,:],dataMatrix[j,:],kernelOption)
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
**
svmQuicken.py
**
import numpy as np
import matplotlib.pyplot as plt
'''
Created on 2017年7月11日
@author: fujianfei
'''
class optStruct:
'''
定義common數(shù)據(jù)結(jié)構(gòu)，存儲(chǔ)需要用到的變量：
'''
def __init__(self, dataMatIn, classLabels, C, toler, kernelOption):
'''
X:訓(xùn)練的數(shù)據(jù)集
labelMat:X對(duì)應(yīng)的類別標(biāo)簽
C:松弛變量系數(shù)
tol:容錯(cuò)率
m:樣本的個(gè)數(shù)
alphas:拉格朗日系數(shù)漓雅，需要優(yōu)化項(xiàng)
b:閾值
eCache:第一列 標(biāo)志位,標(biāo)志Ek是否有效录别，1為有效，0為無效 第二列 錯(cuò)誤率Ek
K:核矩陣
kernelOption:核選項(xiàng)邻吞，如果是線性核kernelOption=('linear', 0) 如果是高斯核kernelOption=('rbf', sigma)组题，sigma為高斯核參數(shù)
'''
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.kernelOpt = kernelOption
self.m,self.n = np.shape(dataMatIn)
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2)))
self.K = np.mat(np.zeros((self.m,self.m)))
#事先把核矩陣都計(jì)算并存儲(chǔ)好，避免以后多次計(jì)算
for i in range(self.m):
self.K[:,i] = calcKernelValue(self.X, self.X[i,:], kernelOption)
def calcEK(oS,k):
'''
計(jì)算誤差Ek
'''
fXk = float(np.multiply(oS.alphas,oS.labelMat).T.dot(oS.K[:,k])) + oS.b
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJ(i,oS,Ei):
'''
啟發(fā)式算法選擇j抱冷，選擇具有最大步長的j
'''
#1.定義步長maxDeltaE (Ei-Ek) 取得最大步長時(shí)的K值maxK 需要返回的Ej （具有最大步長 崔列，即|Ei-Ej|值最大）
maxK = -1; maxDeltaE = 0; Ej = 0
#2.將Ei保存到數(shù)據(jù)結(jié)構(gòu)的eCache中去
oS.eCache[i] = [1,Ei]
#3.定義list validEcacheList，存放有效的Ek
validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
#4.判斷 如果len(validEcacheList)>1 遍歷validEcacheList旺遮，找到最大的|Ei-Ej|
if (len(validEcacheList) > 1):
for k in validEcacheList:
Ek = calcEK(oS, k)
deltaE = abs(Ei - Ek)
if (maxDeltaE < deltaE):
maxDeltaE = deltaE
maxK = k
Ej = Ek
return maxK,Ej
#5.否則就隨機(jī)選擇j
else:
print("---------------隨機(jī)選擇的j---------------------")
j = selectJrand(i,oS.m)
Ej = calcEK(oS, j)
return j,Ej
def updateEk(oS,k):
'''
計(jì)算并更新Ek值到緩存eCache中
'''
Ek = calcEK(oS, k)
oS.eCache[k] = [1,Ek]
def calcfXk(oS,k):
'''
計(jì)算誤差fXk,數(shù)據(jù)集訓(xùn)練結(jié)束后赵讯，可用它來對(duì)testdate進(jìn)行分類
'''
fXk = float(np.multiply(oS.alphas,oS.labelMat).T.dot(oS.K[:,k])) + oS.b
return fXk
def calcKernelValue(matrix_x, sample_x, kernelOption):
'''
定義核函數(shù)
kernelOption=linear 線性
kernelOption=rbf 高斯核函數(shù)
'''
kernelType = kernelOption[0]
numSamples = matrix_x.shape[0]
kernelValue = np.mat(np.zeros((numSamples, 1)))
if kernelType == 'linear':
kernelValue = matrix_x.dot(sample_x.T)
elif kernelType == 'rbf':
sigma = kernelOption[1]
if sigma == 0:
sigma = 1.0
for i in range(numSamples):
diff = matrix_x[i, :] - sample_x
kernelValue[i] = np.exp(diff.dot(diff.T) / (-2.0 * sigma**2))
else:
raise NameError('Not support kernel type! You can use linear or rbf!')
return kernelValue
def selectJrand(i,m):
'''
根據(jù)i，隨機(jī)選擇j
'''
j=i
while (j==i):
j = int(np.random.uniform(0,m))#在0-m的隨機(jī)選一個(gè)數(shù)
return j
def clipAlpha(aj,H,L):
'''
用于調(diào)整大于H或小于L的值耿眉，剪輯最優(yōu)解
'''
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def innerL(i,oS):
'''
內(nèi)循環(huán)边翼，選定i后，在此函數(shù)根據(jù)啟發(fā)式算法選定j,優(yōu)化alphas[i],alphas[j]
計(jì)算優(yōu)化后的Ei,Ej,b跷敬，最后再將它們?nèi)看嫒霐?shù)據(jù)結(jié)構(gòu)optStruct
'''
Ei = calcEK(oS, i)
#判斷優(yōu)化前的alphas[i]是否滿足KKT條件讯私，如果不滿足，進(jìn)行優(yōu)化（啟發(fā)式算法選擇i）
#看論壇上有人文西傀，KKT條件有三個(gè)：alphas[i]=0；alphas[i]=C;0<alphas[i]<C桶癣；而這里只加了0<alphas[i]<C的判斷是不是漏了等于0和等于C的情況
#其實(shí)alphas[i]=0和alphas[i]=C已經(jīng)包含進(jìn)了這個(gè)判斷
#alphas[i]=0時(shí)滿足oS.alphas[i] < oS.C拥褂，故而必要要滿足oS.labelMat[i]*Ei < 0，兩者一和起來不就是alphas[i]=0的KKT條件嗎牙寞。同理饺鹃，alphas[i]=C也是
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ù)啟發(fā)式算法選定j,并計(jì)算好對(duì)應(yīng)的Ej
j,Ej = selectJ(i, oS, Ei)
print("啟發(fā)式算法選出的 i = %d,j= %d" % (i,j))
#重新開辟兩處內(nèi)存，復(fù)制優(yōu)化前的alphas[i]和alphas[j]
#因?yàn)楹竺媾袛鄡?yōu)化后的alphas[j]是否有足夠的變化间雀，需要用到優(yōu)化前的
alphaIold = oS.alphas[i].copy();alphaJold = oS.alphas[j].copy();
#計(jì)算alphas[j]的邊界L,H
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])
#如果最小值L等于最大值H悔详，則沒必要再進(jìn)行優(yōu)化了，直接返回0
if(L == H):print('L=%d == H=%d' % (L,H));return 0
#計(jì)算eta
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
#如過eta>=0惹挟，則可證明最優(yōu)值在邊界處取得茄螃，不需要再優(yōu)化，直接返回0
#解釋：-eta是我們構(gòu)造的拉格朗日函數(shù)L的二階導(dǎo)數(shù)连锯，如果eta>0归苍，二階導(dǎo)數(shù)<0用狱，L在區(qū)間內(nèi)為單調(diào)函數(shù)，所以最優(yōu)值在邊界處取得
#最優(yōu)值解alphas[j]就等于L或H拼弃，此時(shí)不需要優(yōu)化夏伊，也可以將此時(shí)的alphas[j]和對(duì)應(yīng)的alphas[i]保存到oS里
#但一般情況，不會(huì)出現(xiàn)這種情況吻氧，比如我們的線性核函數(shù)溺忧，可以想象成（x1+x2）^2,拆開后2*x1*x2 - x1^2 -X2^2 肯定是<0的。=0的情況太復(fù)雜盯孙，但基本不會(huì)出現(xiàn)鲁森，不考慮。
if(eta >= 0):
print('eta >= 0#################################################################################')
#不需要再優(yōu)化镀梭，直接返回0
return 0
#更新alphas[j]
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
#更新Ej
updateEk(oS, j)
#如果alphas[j]的變化很小刀森，可忽略，則不需再優(yōu)化报账，直接返回0
if (abs(oS.alphas[j]-alphaJold) < 0.00001) :
print('j not moving enough')
#j沒有變化足夠的多研底，不需要再優(yōu)化，直接返回0
return 0
#根據(jù)alphas[j]計(jì)算alphas[i]
oS.alphas[i] += oS.labelMat[i]*oS.labelMat[j]*(alphaJold - oS.alphas[j])
#更新Ei
updateEk(oS, i)
#計(jì)算閾值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[j,i]
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] 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:print("alphas[i]在容錯(cuò)范圍內(nèi)透罢，不需優(yōu)化");return 0
def isFitKKT(oS):
'''
判斷是否在精度范圍內(nèi)符合KKT條件榜晦，符合返回True.作為停機(jī)的最后驗(yàn)證條件.精度為oS.tol
如果符合KKT條件，那么找出的alphas羽圃，一定為最優(yōu)解
'''
for i in range(oS.m):
#如果alphas小于0 或 大于C乾胶，不滿足KKT條件，直接返回False
if (oS.alphas[i] < 0 or oS.alphas[i] > oS.C) : return False
#如果不滿足KKT的核心條件朽寞，之間返回False
if ((oS.alphas[i] == 0 and calcfXk(oS,i) * oS.labelMat[i] < 1) or (oS.alphas[i] > 0 and oS.alphas[i] < oS.C and abs(calcfXk(oS,i) * oS.labelMat[i] - 1) > oS.tol) or (oS.alphas[i] == oS.C and calcfXk(oS,i) * oS.labelMat[i] > 1)) : return False
#如果上面兩個(gè)條件都滿足了识窿，最后再滿足alphas*labelMat之和等于0，則便返回True脑融，符合KKT條件
return abs(oS.alphas.T.dot(oS.labelMat)) < oS.tol
def smoP(dataMatIn,classLabels,C,toler,maxIter,kernelOption):
'''
smo優(yōu)化后的算法
dataMatIn:訓(xùn)練的數(shù)據(jù)集
classLabels:類別標(biāo)簽
C:松弛變量系數(shù)
toler:容錯(cuò)率
kernelOption:核選項(xiàng)喻频，如果是線性核kernelOption=('linear', 0) 如果是高斯核kernelOption=('rbf', sigma)，sigma為高斯核參數(shù)
'''
oS = optStruct(np.mat(dataMatIn),np.mat(classLabels).T,C,toler, kernelOption)#初始化oS
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, paris changes %d" % (iter,i,alphaPairsChanged))
print(isFitKKT(oS))
iter += 1
else :#遍歷非邊界上的值（支持向量機(jī)）
nonBoundIs = np.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, paris changes %d" % (iter,i,alphaPairsChanged))
print(isFitKKT(oS))
iter += 1
if entireSet : entireSet = False
elif (alphaPairsChanged == 0) : entireSet = True
print("iteration number : %d" % iter)
print("-#-#-#-#-#-#-#-#-#_#-#-#_#-#_##_#_#_#_#_#_#_#__#_#_#_#_#_")
print(isFitKKT(oS))
return oS
def showResult(oS):
'''
畫圖
'''
w=np.multiply(oS.alphas,oS.labelMat).T.dot(oS.X)
w=np.mat(w)
x1=oS.X[:,0]
y1=oS.X[:,1]
x2=range(20,100)
b=float(oS.b)
w0=float(w[0,0])
w1=float(w[0,1])
# y2 = [-b/w1-w0/w1*elem for elem in x2]
# plt.plot(x2, y2)
for i in range(oS.m):
if ((oS.alphas[i] < oS.C) and (oS.alphas[i] > 0)):
print('########################')
print(oS.alphas[i])
plt.scatter(x1[i], y1[i],s=60,c='red',marker='o',alpha=0.5,label='SV')
if int(oS.labelMat[i]) == -1:
plt.scatter(x1[i], y1[i],s=30,c='red',marker='.',alpha=0.5,label='-1')
elif int(oS.labelMat[i]) == 1:
plt.scatter(x1[i], y1[i],s=30,c='blue',marker='x',alpha=0.5,label='+1')
plt.show()
def testSVM(svm, test_x, test_y):
'''
測(cè)試訓(xùn)練后結(jié)果正確率
'''
test_x = np.mat(test_x)
test_y = np.mat(test_y).T
numTestSamples = test_x.shape[0]
supportVectorsIndex = np.nonzero(svm.alphas.A > 0)[0]
supportVectors = svm.X[supportVectorsIndex]
supportVectorLabels = svm.labelMat[supportVectorsIndex]
supportVectorAlphas = svm.alphas[supportVectorsIndex]
matchCount = 0
for i in range(numTestSamples):
kernelValue = calcKernelValue(supportVectors, test_x[i, :], svm.kernelOpt)
predict = kernelValue.T.dot(np.multiply(supportVectorLabels, supportVectorAlphas)) + svm.b
if np.sign(predict) == np.sign(test_y[i]):
matchCount += 1
accuracy = float(matchCount) / numTestSamples
return accuracy
init.py
from SVM import svmMLiA,svmQuicken
import numpy as np
if __name__ == '__main__':
dataMatIn,classLabels = svmMLiA.loadDataSet('data2.txt')
C=0.6
toler=0.001
maxIter = 40
#用啟發(fā)式算法版本肘迎，速度快甥温，但是效果不好
oS = svmQuicken.smoP(dataMatIn, classLabels, C, toler, maxIter, ('linear', 0))
#用啟發(fā)式算法版本，速度快妓布，但是效果不好
#不用啟發(fā)式算法姻蚓，速度慢，但是效果好
# b,alphas = svmMLiA.smoSimple(dataMatIn, classLabels, C, toler, maxIter, ('rbf', 0.05))
# oS = svmQuicken.optStruct(np.mat(dataMatIn),np.mat(classLabels).T,C,toler, ('rbf', 0.05))#初始化oS
# oS.alphas = alphas
# oS.b = b
#不用啟發(fā)式算法匣沼，速度慢狰挡，但是效果好
print(oS.alphas)
#計(jì)算正確率
rightRate = svmQuicken.testSVM(oS,dataMatIn, classLabels)
#畫圖
svmQuicken.showResult(oS)
print("正確率是 : %f" % rightRate)
參考文獻(xiàn)：
[1]李航的《統(tǒng)計(jì)學(xué)習(xí)方法》，清華大學(xué)出版社
