实验背景
SVM算法在上世纪60年代就已经被提出，学名为Support Vector Machine，是一种非常经典的监督学习方法。我在看来，SVM是最好的现成的分类器，这里说的“现成”指的是分类器不加修改即可直接使用。这意味着在数据上应用基本形式（没有针对数据进行修改）的SVM分类器就可以得到低错误率的结果。SVM能够对训练集之外的数据点做出很好的分类决策。
SVM有很多实现方法，本篇博客只介绍最流行的一种方法，序列最小优化（Sequential Minimal Optimization,SMO）算法。
1.SVM算法原理
1.1.分隔超平面
首先我们要介绍下SVM支持向量机的目的，如图所示
对于整个图像来说，如何找到一条直线将+，-，完美的分开是很困难的，而这条将数据集分开的直线称为分隔超平面。图像中的数据点都在二维平面上，所以分隔超平面就是一条直线，如果数据集是三维的，那分隔超平面就是一个平面。更高维的情况以此类推。SVM的目的就是找到这个最佳的分隔超平面。而SVM中的支持向量，就是指离分隔超平面最近的那些点，而这些点离分隔面越大的话，这个分隔面就越接近完美分隔面。
1.2.最大间隔
1.2.1.寻找最小间隔数据点
如图所示，对于超平面来说，距离其最近的点就是被圆圈框起来的一个+号和两个-号。
假设这个超平面方程为：
那么位于超平面上方的支持向量所经过的直线可以表达为：
位于超平面下方的支持向量所经过的直线可以表达为：
根据直线之间的距离计算公式，最后可以化为：
将问题化为寻找参数w和b，使得下述公式最大
1.2.2.拉格朗日乘子法
为了更好的求得最大值，我们引入拉格朗日乘子得到对应的拉格朗日函数
然后，令L（w,b,α）对w和b的偏导为零
将w，b带回原式得
等价形式为
最终模型：
其中约束条件为：
1.3.SMO高效优化算法
SMO算法的工作原理是：每次循环中选择两个合适的α进行优化处理。一旦找到一对合适的α，就增大其中一个的同时减小另一个。
流程图大致如下：
公式表达为：
得到最优解为：
带回原式解得w和b为
得到分类平面：
伪代码大致如下：
创建一个α向量并将其初始化为0向量
当迭代次数小于最大迭代次数时（外循环）
对数据集中的每个数据向量（内循环）；
如果该数据向量可以被优化：
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化，退出内循环
如果所有向量都没被优化，增加迭代数目，继续下一次循环
2.实验结果
from time import sleep
import matplotlib.pyplot as plt
import numpy as np
import random
import types
"""
函数说明:读取数据
输入：
fileName - 文件名
输出：
dataMat - 数据矩阵
labelMat - 数据标签
"""
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
"""
函数说明:数据可视化
输入：
dataMat - 数据矩阵
labelMat - 数据标签
输出：
无
"""
def showDataSet(dataMat, labelMat):
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) #转换为numpy矩阵
data_minus_np = np.array(data_minus) #转换为numpy矩阵
plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1]) #正样本散点图
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1]) #负样本散点图
plt.show()
"""
函数说明:随机选择alpha
输入：
i - alpha
m - alpha参数个数
输出：
j -
"""
def selectJrand(i, m):
j = i #选择一个不等于i的j
while (j == i):
j = int(random.uniform(0, m))
return j
"""
函数说明:修剪alpha
输入：
aj - alpha值
H - alpha上限
L - alpha下限
输出：
aj - alpah值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
函数说明:简化版SMO算法
输入：
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
输出：
无
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
#转换为numpy的mat存储
dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
#初始化b参数，统计dataMatrix的维度
b = 0; m,n = np.shape(dataMatrix)
#初始化alpha参数，设为0
alphas = np.mat(np.zeros((m,1)))
#初始化迭代次数
iter_num = 0
#最多迭代matIter次
while (iter_num < maxIter):
alphaPairsChanged = 0
for i in range(m):
#步骤1：计算误差Ei
fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
Ei = fXi - float(labelMat[i])
#优化alpha，更设定一定的容错率。
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
#随机选择另一个与alpha_i成对优化的alpha_j
j = selectJrand(i,m)
#步骤1：计算误差Ej
fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j])
#保存更新前的aplpha值，使用深拷贝
alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
#步骤2：计算上下界L和H
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
#步骤3：计算eta
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
#步骤4：更新alpha_j
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
#步骤5：修剪alpha_j
alphas[j] = clipAlpha(alphas[j],H,L)
if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
#步骤6：更新alpha_i
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
#步骤7：更新b_1和b_2
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
#步骤8：根据b_1和b_2更新b
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
class optStruct:
"""
数据结构，维护所有需要操作的值
Parameters：
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
kTup - 包含核函数信息的元组,第一个参数存放核函数类别，第二个参数存放必要的核函数需要用到的参数
"""
def __init__(self, dataMatIn, classLabels, C, toler, kTup):
self.X = dataMatIn #数据矩阵
self.labelMat = classLabels #数据标签
self.C = C #松弛变量
self.tol = toler #容错率
self.m = np.shape(dataMatIn)[0] #数据矩阵行数
self.alphas = np.mat(np.zeros((self.m,1))) #根据矩阵行数初始化alpha参数为0
self.b = 0 #初始化b参数为0
self.eCache = np.mat(np.zeros((self.m,2))) #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。
self.K = np.mat(np.zeros((self.m,self.m))) #初始化核K
for i in range(self.m): #计算所有数据的核K
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
def kernelTrans(X, A, kTup):
"""
通过核函数将数据转换更高维的空间
Parameters：
X - 数据矩阵
A - 单个数据的向量
kTup - 包含核函数信息的元组
输出：
K - 计算的核K
"""
m,n = np.shape(X)
K = np.mat(np.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 = np.exp(K/(-1*kTup[1]**2)) #计算高斯核K
else: raise NameError('核函数无法识别')
return K #返回计算的核K
def calcEk(oS, k):
"""
计算误差
Parameters：
oS - 数据结构
k - 标号为k的数据
输出：
Ek - 标号为k的数据误差
"""
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJ(i, oS, Ei):
"""
内循环启发方式2
Parameters：
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
输出：
j, maxK - 标号为j或maxK的数据的索引值
Ej - 标号为j的数据误差
"""
maxK = -1; maxDeltaE = 0; Ej = 0 #初始化
oS.eCache[i] = [1,Ei] #根据Ei更新误差缓存
validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] #返回误差不为0的数据的索引值
if (len(validEcacheList)) > 1: #有不为0的误差
for k in validEcacheList: #遍历,找到最大的Ek
if k == i: continue #不计算i,浪费时间
Ek = calcEk(oS, k) #计算Ek
deltaE = abs(Ei - Ek) #计算|Ei-Ek|
if (deltaE > maxDeltaE): #找到maxDeltaE
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej #返回maxK,Ej
else: #没有不为0的误差
j = selectJrand(i, oS.m) #随机选择alpha_j的索引值
Ej = calcEk(oS, j) #计算Ej
return j, Ej #j,Ej
def updateEk(oS, k):
"""
计算Ek,并更新误差缓存
输入：
oS - 数据结构
k - 标号为k的数据的索引值
输出：
无
"""
Ek = calcEk(oS, k) #计算Ek
oS.eCache[k] = [1,Ek] #更新误差缓存
def innerL(i, oS):
"""
优化的SMO算法
输入：
i - 标号为i的数据的索引值
oS - 数据结构
输出：
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
#步骤1：计算误差Ei
Ei = calcEk(oS, i)
#优化alpha,设定一定的容错率。
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)):
#使用内循环启发方式2选择alpha_j,并计算Ej
j,Ej = selectJ(i, oS, Ei)
#保存更新前的aplpha值，使用深拷贝
alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
#步骤2：计算上下界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])
if L == H:
print("L==H")
return 0
#步骤3：计算eta
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
if eta >= 0:
print("eta>=0")
return 0
#步骤4：更新alpha_j
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
#步骤5：修剪alpha_j
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
#更新Ej至误差缓存
updateEk(oS, j)
if (abs(oS.alphas[j] - alphaJold) < 0.00001):
print("alpha_j变化太小")
return 0
#步骤6：更新alpha_i
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
#更新Ei至误差缓存
updateEk(oS, i)
#步骤7：更新b_1和b_2
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]
#步骤8：根据b_1和b_2更新b
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
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin',0)):
"""
输入：
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
kTup - 包含核函数信息的元组
输出：
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler, kTup) #初始化数据结构
iter = 0 #初始化当前迭代次数
entireSet = True; alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): #遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环
alphaPairsChanged = 0
if entireSet: #遍历整个数据集
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS) #使用优化的SMO算法
print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
else: #遍历非边界值
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] #遍历不在边界0和C的alpha
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet: #遍历一次后改为非边界遍历
entireSet = False
elif (alphaPairsChanged == 0): #如果alpha没有更新,计算全样本遍历
entireSet = True
print("迭代次数: %d" % iter)
return oS.b,oS.alphas #返回SMO算法计算的b和alphas
def img2vector(filename):
returnVect = np.zeros((1, 1024))
fr = open(filename)
for i in range(32):
lineStr = fr.readline()
for j in range(32):
returnVect[0, 32 * i + j] = int(lineStr[j])
return returnVect
def loadImages(dirName):
from os import listdir
hwLabels = []
trainingFileList = listdir(dirName)
m = len(trainingFileList)
trainingMat = np.zeros((m, 1024))
for i in range(m):
fileNameStr = trainingFileList[i]
fileStr = fileNameStr.split('.')[0]
classNumStr = int(fileStr.split('_')[0])
if classNumStr == 9:
hwLabels.append(-1)
else:
hwLabels.append(1)
trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
return trainingMat, hwLabels
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadImages('D:/vscode/python/.vscode/trainingDigits')
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
datMat = np.mat(dataArr);
labelMat = np.mat(labelArr).transpose()
svInd = np.nonzero(alphas.A > 0)[0]
sVs = datMat[svInd]
labelSV = labelMat[svInd];
print("支持向量个数： %d " % np.shape(sVs)[0])
m, n = np.shape(datMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("训练集错误率: %f" % (float(errorCount) / m))
dataArr, labelArr = loadImages('D:/vscode/python/.vscode/testDigits')
errorCount = 0
datMat = np.mat(dataArr);
labelMat = np.mat(labelArr).transpose()
m, n = np.shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("测试集错误率: %f" % (float(errorCount) / m))
if __name__ == '__main__':
testDigits()
支持向量数量有点过多，应该可以进一步优化
3.总结
1.SVM支持向量机是一种分类器，它的泛化错误率低，具有良好的学习能力，且得到的结果具有很好的推广性。
2.SVM支持向量机通过求解一个二次优化问题来最大化分类间隔。相比最初始版本采用的复杂且低效的二次规划求解方法，SMO算法通过优化求解两个α值来提高训练速度。同时还有进一步提高的空间。