支持向量機的理論支持在此不細說，可以參考李航的《統計學習》，還有西瓜書。

簡化版SMO算法處理小規模數據集

SMO算法是一種啟發式算法。此簡化版首先在數據集上遍歷每一個alpha，然后在剩下的alpha集合中隨機選擇另一個alpha，從而建立alpha對。

# -*- coding: utf-8 -*- from numpy import * from time import sleep# SMO算法的輔助函數 def loadDataSet(fileName): #加載并預處理數據集dataMat = []; labelMat = []fr = open(fileName,'r')for line in fr.readlines():lineArr = line.strip().split('\t') # 以制表符分割dataMat.append([float(lineArr[0]), float(lineArr[1])]) #提取前兩個元素存入data.Mat中labelMat.append(float(lineArr[2])) # [].append(),最終的形式是矩陣return dataMat,labelMatdef selectJrand(i,m): # 該輔助函數用于在某個區間范圍內隨機選擇一個整數j=i # m是所有alpha的數目，i是第一個alpha的下標while (j==i):j = int(random.uniform(0,m)) # random.uniform(0,m)用于生成指定范圍內的隨機浮點數return jdef clipAlpha(aj,H,L): # 該輔助函數用于在數值太大時對其進行調整if aj > H: aj = Hif L > aj:aj = Lreturn aj# 簡化版SMO算法 def smoSimple(dataMatIn, classLabels, C, toler, maxIter): # 參數：數據集，類別標簽，常數c，容錯率，循環次數dataMatrix = mat(dataMatIn) # mat()轉換成矩陣類型labelMat = mat(classLabels).transpose() #轉置之前是列表，轉置后是一個列向量b = 0; m,n = shape(dataMatrix) # 得到行，列數，m行，n列alphas = mat(zeros((m,1))) # zeros(shape, dtype=float, order='C'),所以也可以寫作zeros((10,1),)iter = 0 # 該變量存儲的是在沒有任何alpha改變時遍歷數據集的次數while (iter < maxIter): # 限制循環迭代次數，也就是在數據集上遍歷maxIter次，且不再發生任何alpha修改，則循環停止alphaPairsChanged = 0 # 每次循環時先設為0，然后再對整個集合順序遍歷，該變量用于記錄alpha是否已經進行優化for i in range(m): # 遍歷每行數據向量，m行# 該公式是分離超平面，我們預測值fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b #print 'fxi:',fxiEi = fXi - float(labelMat[i]) # 預測值和真實輸出之差 # 如果誤差很大就對該數據對應的alpha進行優化，正負間隔都會被測試，同時檢查alpha值 if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): j = selectJrand(i,m) # 隨機選擇不等于i的0-m的第二個alpha值fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + bEj = fXj - float(labelMat[j])alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();if (labelMat[i] != labelMat[j]): # 這里是對SMO最優化問題的子問題的約束條件的分析L = max(0, alphas[j] - alphas[i]) # L和H分別是alpha所在的對角線端點的界H = min(C, C + alphas[j] - alphas[i]) # 調整alphas[j]位于0到c之間else:L = max(0, alphas[j] + alphas[i] - C)H = min(C, alphas[j] + alphas[i])if L==H: print "L==H"; continue # L=H停止本次循環# 是一個中間變量：eta=2xi*xi-xixi-xjxj，是alphas[j]的最優修改量eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T \- dataMatrix[j,:]*dataMatrix[j,:].Tif eta >= 0: print "eta>=0"; continue # eta>=0停止本次循環，這里是簡化計算alphas[j] -= labelMat[j]*(Ei - Ej)/eta # 沿著約束方向未考慮不等式約束時的alpha[j]的解alphas[j] = clipAlpha(alphas[j],H,L) # 此處是考慮不等式約束的alpha[j]解if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue # 如果該alpha值不再變化，就停止該alpha的優化alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) # 更新alpha[i]# 完成兩個alpha變量的更新后，都要重新計算閾值bb1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T \- labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]* dataMatrix[j,:].T #李航統計學習7.115式b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T \- labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T #李航統計學習7.116式if (0 < alphas[i]) and (C > alphas[i]): b = b1elif (0 < alphas[j]) and (C > alphas[j]): b = b2 else: b = (b1 + b2)/2.0 # alpha[i]和alpha[j]是0或者c,就取中點作為balphaPairsChanged += 1 # 到此的話說明已經成功改變了一對alphaprint "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)if (alphaPairsChanged == 0): iter += 1 # 如果alpha不再改變迭代次數就加1else: iter = 0 print "iteration number: %d" % iterreturn b,alphas# 主函數 dataArr,labelArr=loadDataSet('testSet.txt') # 因為在同一文件夾下，就不用寫絕對路徑 #print dataArr b,alphas=smoSimple(dataArr, labelArr, 0.6, 0.001, 40) print 'b:',b print 'alphas[alphas>0]:',alphas[alphas>0] # 數組過濾 print shape(alphas[alphas>0]) # 得到支持向量的個數 for i in range(100): # 得到是支持向量的數據點if alphas[i]>0.0: print dataArr[i],labelArr[i]

其中要注意的python語法：

• 數組和矩陣的轉換

In [17]:type(dataArr) # 轉換之前是list的類型 Out[17]: listIn [18]:mat(dataArr) # 轉換之后是矩陣類型 Out[18]: matrix([[ 3.542485, 1.977398],[ 3.018896, 2.556416],[ 7.55151 , -1.58003 ],..., [ 2.912122, -0.202359],[ 1.731786, 0.589096],[ 2.387003, 1.573131]])

• Python自帶的copy()，deepcopy()，numpy的copy（）之間的區別

>>> import copy >>> origin = [1, 2, [3, 4]] >>> cop1 = copy.copy(origin) # 淺拷貝 >>> cop2 = copy.deepcopy(origin) # 深拷貝 >>> cop1 == cop2 True >>> cop1 is cop2 False #cop1 和 cop2 看上去相同，但已不再是同一個object >>> origin[2][0] = "hey!" >>> origin [1, 2, ['hey!', 4]] >>> cop1 [1, 2, ['hey!', 4]] >>> cop2 [1, 2, [3, 4]] #把origin內的子list [3, 4] 改掉了一個元素，觀察 cop1 和 cop2

可以看到 cop1，也就是 shallow copy 跟著 origin 改變了。而 cop2 ，也就是 deep copy 并沒有變。
似乎 deep copy 更加符合我們對「復制」的直覺定義: 一旦復制出來了，就應該是獨立的了。如果我們想要的是一個字面意義的「copy」，那就直接用 deep_copy 即可。

注意：這里指的是python自帶的copy( )函數，copy.copy(object)是淺拷貝，而numpy的copy( )函數，即object.copy( )是深拷貝的效果。

from numpy import * In [26]:origin Out[26]: matrix([[1, 2], [3, 4]])In [28]:origin[0][0,1] Out[28]: 2In [29]:origin[0][0,0] Out[29]: 1 In [31]:old=origin[0][0,1].copy()In [32]:origin[0][0,1]=5In [33]:origin[0][0,1] Out[33]: 5In [34]:old Out[34]: 2

分類圖示：

在園圈中的就是支持向量

代碼中的公式含義：

（1）

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

fxi對應的是《統計學習》中7.104的g(x)：

g(x)=∑i=1Nα?iyixi?x+b?
其代表的是分離超平面
（2）

Ei = fXi - float(labelMat[i])

Ei=g(xi)?yi
表示預測值和真實輸出之差
（3）

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])

L和H是alphas的所在對角線端點的界：
如果y1≠y2則L=max(0,αoldj?αoldi),H=min(c,c+αoldj?αoldi)
如果y1=y2則L=max(0,αoldj+αoldi?c),H=min(c,αoldj+αoldi)
（4）

eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T \- dataMatrix[j,:]*dataMatrix[j,:].T

對應

η=2xj?xi?xi?xi?xj?xj
表示的是alphas[j]的最優修改量，是一個中間變量
（5）

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

對應

αnew,uncj=αoldj?yj(Ei?Ej)η
是沿著約束方向未考慮不等式約束時的alpha[j]的解

（6）

alphas[j] = clipAlpha(alphas[j],H,L)

此處是考慮不等式約束的alpha[j]解：

αnewj=???????Hαnew,uncjLαnew,uncj>HL?αnew,uncj?Hαnew,uncj<L
得到了第一個變量的值
（7）

alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])

αnewi=αoldi+y1y2(αoldj?αnewj)

（8）

b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T \- labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]* dataMatrix[j,:].T

得到的是b：

bnew1=bold?Ei?yixi?xi(αnewi?αoldi)?y2xj?xi(αnewj?αoldj)+bold

代碼中的主要公式大概就這些。

運行結果：

j not moving enough j not moving enough iteration number: 1 j not moving enough j not moving enough ..., j not moving enough iteration number: 38 j not moving enough j not moving enough iteration number: 39 j not moving enough j not moving enough iteration number: 40 b: [[-3.83495394]] alphas[alphas>0]: [[ 0.15601002 0.14181599 0.06893826 0.36676427]] (1L, 3L) [4.658191, 3.507396] -1.0 [3.457096, -0.082216] -1.0 [6.080573, 0.418886] 1.0

“不忘初心，方得始終”搞了那么多代碼和公式，不要沉浸進去而不知道最后要求的是什么？SMO算法的目標是求出一系列的α和b，而此時運行得到的就是那些支持向量和相對應的參數α，還有b，由此也就得到了分離超平面，就是：

∑i=1Nα?iyixi?x+b?=0
分類決策函數可以寫成：
f(x)=sign(∑i=1Nα?iyixi?x+b?)
由此驗證了在決定分離超平面時只有支持向量起作用。如果移動支持向量將改變所求的解；但如果在間隔邊界以外移動其他實例點，甚至去掉這些點，則解不會改變。
由于SMO算法的隨機性，每次的運行結果可能不同。

還有一點要注意：就是優化結束的同時必須確保合適的時機結束循環，如果程序執行到for循環的最后一行都不執行continue語句，那么就成功地改變了一對α值，同時可以增加alphaPairsChanged的值，在for循環之外，需要檢查α值是否有更新，如果有更新則將iter設為0后繼續運行程序。只有在所有的數據集上遍歷maxIter次，且不再發生任何α的修改后，程序才停止并退出while循環。

利用完整版platt SMO算法加速優化

在完整版的SMO算法中，實現alpha的更改和代數運算的優化環節一模一樣，在優化過程中，唯一不同的是選擇alpha的方式，完整版的應用了一些能夠提速的啟發式方法。

# -*- coding: utf-8 -*- """ Created on Wed Oct 18 20:53:40 2017@author: LiLong """ from numpy import * from time import sleep# SMO算法的輔助函數 def loadDataSet(fileName): #加載并預處理數據集dataMat = []; labelMat = []fr = open(fileName,'r')for line in fr.readlines():lineArr = line.strip().split('\t') # 以制表符分割dataMat.append([float(lineArr[0]), float(lineArr[1])]) #提取前兩個元素存入data.Mat中labelMat.append(float(lineArr[2])) # [].append(),最終的形式是矩陣return dataMat,labelMat# 完整版platt SMO 算法的支持函數 # 建立一個數據結構來保存所有的重要值 class optStruct:def __init__(self,dataMatIn, classLabels, C, toler): self.X = dataMatInself.labelMat = classLabelsself.C = Cself.tol = tolerself.m = shape(dataMatIn)[0] # 有多少行數據self.alphas = mat(zeros((self.m,1)))self.b = 0self.eCache = mat(zeros((self.m,2))) # 誤差緩存，第一列是ecache是否有效的標志位，第二列是實際的E值def clipAlpha(aj,H,L):if aj > H: aj = Hif L > aj:aj = Lreturn ajdef selectJrand(i,m): # 該輔助函數用于在某個區間范圍內隨機選擇一個整數j=i # m是所有alpha的數目，i是第一個alpha的下標while (j==i):j = int(random.uniform(0,m)) # random.uniform(0,m)用于生成指定范圍內的隨機浮點數return j# 計算E值并返回,E值是函數對輸入xi的預測值與真實輸出的差 def calcEk(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T))+ oS.bEk = fXk - float(oS.labelMat[k])return Ek# 用于選擇合適的第二個alpha值以保證每次優化中采用最大步長，是內循環的啟發式方法 def selectJ(i, oS, Ei): # 該函數的誤差值與第一個alpha值Ei和下標i有關maxK = -1; maxDeltaE = 0; Ej = 0oS.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) # 傳遞對象和k，計算誤差值deltaE = abs(Ei - Ek)if (deltaE > maxDeltaE): # 選擇具有最大步長的jmaxK = k; maxDeltaE = deltaE; Ej = Ek # 會在所有的值上循環，并選擇其中使得改變最大的那個值return maxK, Ejelse: # 在這種情況下（第一次，我們沒有任何有效的eCache值 ），隨機選擇一個alpha值j = selectJrand(i, oS.m) Ej = calcEk(oS, j)return j, Ejdef updateEk(oS, k): # alpha改變時更新緩存中的值Ek = calcEk(oS, k) oS.eCache[k] = [1,Ek]# 完整platt 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] > 0)):j,Ej = selectJ(i, oS, Ei) # 第二個alpha選擇中的啟發式方法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 0oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/etaoS.alphas[j] = clipAlpha(oS.alphas[j],H,L)updateEk(oS, j) # 更新誤差緩存if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0oS.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,:].Tb2 = 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,:].Tif (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2else: oS.b = (b1 + b2)/2.0return 1 # 如果有任意一對alpha發生改變，那么就會返回1，其他返回0else: return 0 # 完整版platt SMO的外循環代碼 def smoP(dataMatIn, classLabels, C, toler, maxIter): # 建立一個數據結構來容納所有的數據oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler) iter = 0entireSet = True; alphaPairsChanged = 0 # 退出循環的變量的一些初始化# 迭代次數超過指定的最大值或者遍歷整個集合都未對任意的alpha對進行修改時就退出循環while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): alphaPairsChanged = 0 if entireSet: for i in range(oS.m): # 一開始在數據集上遍歷任意可能的alpha # 選擇第二個alpha，并在可能時對其進行優化處理,有任一一對alpha發生變化化了alphaPairsChanged+1alphaPairsChanged += innerL(i,oS)print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)iter += 1else: # 遍歷所有的非邊界alpha值，也就是不在邊界0或c上的值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 += 1if entireSet: entireSet = False # 在非邊界循環和完整遍歷之間進行切換elif (alphaPairsChanged == 0): entireSet = True print "iteration number: %d" % iterreturn oS.b,oS.alphas # 分類超平面的w計算 def calcWs(alphas,dataArr,classLabels):X = mat(dataArr); labelMat = mat(classLabels).transpose()m,n = shape(X)w = zeros((n,1))for i in range(m):w += multiply(alphas[i]*labelMat[i],X[i,:].T)return w# 主函數 dataArr,labelArr=loadDataSet('testSet.txt') b,alphas=smoP(dataArr,labelArr,0.6,0.001,40) print 'b:',b print 'alphas:',alphas # 輸出w和b ws=calcWs(alphas,dataArr,labelArr) print 'ws:',ws datmat=mat(dataArr) result=datmat[0]*mat(ws)+b # 進行分類 print 'result',result

一些Python的技巧：
nonzero的用法
matrix.A用法

In [5]:a=mat(zeros((3,2))) In [6]:a Out[6]: matrix([[ 0., 0.],[ 0., 0.],[ 0., 0.]]) In [7]:nonzero(a[:,0].A)[0] Out[7]: array([], dtype=int64) In [8]:a[:,0].A Out[8]: array([[ 0.],[ 0.],[ 0.]]) In [15]: a=mat([[0,1],[1,0],[0,1]]) In [16]: a[:,0].A Out[16]: array([[0],[1],[0]]) In [19]: nonzero(a[:,0].A)[0] Out[19]: array([1], dtype=int64) In [20]: a=mat([[0,1],[1,0],[1,1]])In [21]: a[:,0].A Out[21]: array([[0],[1],[1]])In [22]: nonzero(a[:,0].A)[0] Out[22]: array([1, 2], dtype=int64) >>> from numpy import * >>> aa=[0,1,2,3,4,10,6,7,8,0] >>> b=mat(aa) >>> b matrix([[ 0, 1, 2, 3, 4, 10, 6, 7, 8, 0]]) >>> b.A array([[ 0, 1, 2, 3, 4, 10, 6, 7, 8, 0]]) >>> b.A>3 array([[False, False, False, False, True, True, True, True, True,False]], dtype=bool) >>> b.A<4 array([[ True, True, True, True, False, False, False, False, False,True]], dtype=bool) >>> (b.A>3)*(b.A<4) array([[False, False, False, False, False, False, False, False, False,False]], dtype=bool) >>> nonzero((b.A>3)*(b.A<4)) (array([], dtype=int64), array([], dtype=int64)) >>> nonzero((b.A>3)*(b.A<4))[0] array([], dtype=int64) >>>

SMO中拉格朗日乘子的啟發式選擇方法：
所謂的啟發式選擇方法主要思想是每次選擇拉格朗日乘子的時候，優先選擇樣本前面系數0<ai<c作優化（論文中稱為無界樣例），因為在界上（ai為0或C）的樣例對應的系數ai一般不會更改。

這條啟發式搜索方法是選擇第一個拉格朗日乘子用的，那么這樣選擇的話，是否最后會收斂。可幸的是Osuna定理告訴我們只要選擇出來的兩個ai中有一個違背了KKT條件，那么目標函數在一步迭代后值會減小。違背KKT條件不代表0<ai<c=0后，先對所有樣例進行循環，循環中碰到違背KKT條件的（不管界上還是界內）都進行迭代更新。等這輪過后，如果沒有收斂，第二輪就只針對0<ai<c，選擇第二個乘子能夠最大化|E1?E2|。即當E1為正時選擇負的絕對值最大的E2，反之，選擇正值最大的E2。

最后的收斂條件是在界內（0<ai<c只在極小的范圍內變動。

參考：

http://www.cnblogs.com/jerrylead/archive/2011/03/18/1988419.html#3793878

運行結果：

L==H fullSet, iter: 0 i:0, pairs changed 0 L==H fullSet, iter: 0 i:1, pairs changed 0 fullSet, iter: 0 i:2, pairs changed 1 L==H ..., j not moving enough fullSet, iter: 2 i:97, pairs changed 0 fullSet, iter: 2 i:98, pairs changed 0 fullSet, iter: 2 i:99, pairs changed 0 iteration number: 3 b: [[-2.89901748]] alphas: [[ 0.06961952][ 0. ][ 0. ]..., [ 0. ][ 0. ][ 0. ]] ws: [[ 0.65307162][-0.17196128]] result [[-0.92555695]]

得到α和b后就可以進行分類了：
sign(∑i=1Nα?iyixi?x+b?)

完整版的SMO算法和簡化版的幾點區別：

• while的循環退出條件更多一些
• maxiter變量的簡化版的作用不同，簡化版的是當沒有任何alpha發生變化時會將整個集合的一次遍歷過程計成一次迭代，而完整版的的一次迭代計成定義為一次循環過程，而不管該循環具體做了什么事
• while循環內部和簡化版的有所不同
• 速度有很大提升

支持向量機有些復雜，以后在應用中再進一步學習！！

總結

以上是生活随笔為你收集整理的支持向量机—SMO算法源码分析（1）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。