Platt SMO 支持向量机算法（Python实现）

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