支持向量机
优点:泛化错误率低,计算开销不大,结果易解释 缺点:对参数调节和核函数的选择敏感,原始分类器不加修改仅适用于处理二类问题 适用数据类型:数值型和标称型数据
??(有人说)SVM是最好的现成的分类器,“现成”是指分类器不加修饰即可直接使用。这也就意味着数据上应用基本的SVM分类器就可以得到低错误率的结果。SVM对训练集之外的数据点做出很好的分类决策。
基于最大间隔的分隔数据
??两组数据分隔的足够开,很容易用一条直线将两组数据点分开,这组数据被称为线性可分数据,将数据集分隔开来的直线称为分割超平面。因此此时的分割超平面是一条直线。若数据集是三维的,此时用来分隔的是一个平面,更高维度的情况以此类推…n维数据集需要n-1维的对象来进行分割,该对象称为超平面,也就是决策的边界。 ??支持向量就是离分隔超平面最近的那些点,接下来要试着最大化支持向量到分隔面的距离,需要找到此问题的最优化求解方法。
寻找最大间隔
??分割超平面的形式可以写成
w
T
x
+
b
w^Tx+b
wTx+b。计算A点到分隔超平面的距离,必须给出分隔面的法线或垂直的长度,该值为:
∣
w
T
x
+
b
∣
?
1
∥
w
∥
|w^Tx+b|\cdot \dfrac {1}{\parallel w\parallel}
∣wTx+b∣?∥w∥1?,这里的常熟b,类似于Logistics回归中的截距
w
w
w。
1、分类器求解的优化问题
??分类器工作原理:输入数据给分类器会输出一个类别的标签,这相当于是一个类似于Sigmoid的函数的作用。下面使用类似海维赛德阶跃函数(单位阶跃函数)最函数
w
T
x
+
b
w^Tx+b
wTx+b作用得到
f
(
w
T
x
+
b
)
f(w^Tx+b)
f(wTx+b),其中当u<0时
f
(
u
)
f(u)
f(u) 输出-1,反之输出+1。 ??当计算点到分隔面的距离并确定分隔面的位置时,间隔通过
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)来计算,若数据点处于正方向(即+1类)且距离分割超平面很远的位置时,
w
T
x
+
b
w^Tx+b
wTx+b是一个很大的正数,同时
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)也是一个很大的正数。若数据点处于负方向(即-1类)且距离分割超平面很远的位置时,分类标签为-1,
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)也是一个很大的正数。 ??现在目标就是找到
w
w
w和
b
b
b。先找到具有最小间隔的数据点,这些数据点也就是前面提到的支持向量。找到具有最小间隔的数据点需要对该间隔最大化,可以写作:
a
r
g
?
m
a
x
w
,
b
{
m
i
n
n
(
l
a
b
e
l
?
(
w
T
x
+
b
)
)
?
1
∥
w
∥
}
arg\, \mathop{max}\limits_{w,b} \{ \mathop{min}\limits_{n}(label\cdot (w^Tx+b))\cdot \dfrac {1}{\parallel w\parallel} \}
argw,bmax?{nmin?(label?(wTx+b))?∥w∥1?} ??若令所有支持向量的
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)都为1,那么可以通过计算
∥
w
∥
?
1
\parallel w\parallel ^{-1}
∥w∥?1的最大值来得到最终的解。但是并非所有数据点的
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)都为1,只有那些离分隔超平面最近的点得到的值才为1。而离超平面越远的数据点
l
a
b
e
l
?
(
w
T
x
+
b
)
label*(w^Tx+b)
label?(wTx+b)值越大。 ??因此该问题是一个带约束条件的优化问题,约束条件是
l
a
b
e
l
?
(
w
T
x
+
b
)
≥
1.0
label*(w^Tx+b)\geq 1.0
label?(wTx+b)≥1.0。对于这类问题有一个非常著名的求解方法,即拉格朗日乘子法。通过引入拉格朗日乘子,可以基于约束条件来表达原来的问题。由于约束条件都是基于数据点的,因此可以将超平面写成数据点的形式,最终写成:
m
a
x
a
[
∑
i
=
1
m
a
?
1
2
∑
i
,
j
=
1
m
l
a
b
e
l
(
i
)
?
l
a
b
e
l
(
j
)
?
a
i
?
a
j
?
x
(
i
)
,
x
(
j
)
?
]
\mathop{max}\limits_{a}[\sum\limits_{i=1}^m a -\dfrac{1}{2}\sum\limits_{i,j=1}^m label^{(i)}\cdot label^{(j)}\cdot a_i\cdot a_j \langle x^{(i)},x^{(j)} \rangle]
amax?[i=1∑m?a?21?i,j=1∑m?label(i)?label(j)?ai??aj??x(i),x(j)?] 其约束条件是
a
≥
0
,
∑
i
=
1
m
a
i
?
l
a
b
e
l
(
i
)
=
0
a\geq 0 \quad ,\quad \sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
a≥0,i=1∑m?ai??label(i)=0 ??但是这里有一个假设:数据必须线性可分。松弛变量:允许有些数据点可以处于分隔面的错误一侧,这样优化目标能保持不变,此时约束条件变为:
C
≥
a
≥
0
,
∑
i
=
1
m
a
i
?
l
a
b
e
l
(
i
)
=
0
C \geq a\geq 0 \quad ,\quad \sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
C≥a≥0,i=1∑m?ai??label(i)=0 ??这里的常量C用于控制“最大化间隔”和“保证大部分数据点的函数间隔小于1.0”这两个目标的权重。一旦求出了所有的alpha,那么分隔超平面就可以通过这些alpha来表达。SVM的主要工作就是求解这些alpha。
2、SVM应用的一般框架
SVM的一般流程
- 收集数据:可以使用任何方法
- 准备数据:需要数值型数据
- 分析数据:有助于可视化分隔超平面
- 训练算法:SVM的大部分时间都源于训练,该过程主要实现两个参数的调优
- 测试算法:简单的计算过程就可以实现
- 使用算法:几乎所有分类问题都可以使用SVM,SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改
SMO高效优化算法
??下面对两个式子进行优化:一个是最小化的目标函数,另一个是在优化过程中必须遵守的约束条件。
1、Platt的SMO算法
??SMO表示序列最小优化,Platt的SMO算法是将大优化算法分解成小优化问题来求解。这些小优化问题往往很容易求解,并且对他们进行顺序求解的结果与将他们作为整体来求解的结果完全一致。在结果完全一致时,SMO算法求解时间短很多。 ??SMO算法的目标是求出一系列的alpha和b,一旦求出了这些alpha就很容易求出权重向量w并得到分隔超平面。 ??SMO算法的工作原理:每次循环中选择两个alpha进行优化,一旦找到一对合适的alpha就增大其中一个,同时减小另一个。
2、应用简化版SMO算法处理小规模数据集
??Platt SMO算法的完整实现要大量的代码,所以先实现简化版,过程:首先在数据集上遍历所有的alpha,然后在剩下的alpha集合中随机选择另一个alpha,从而构建alpha对。我们要同时改变两个alpha,因为有一个约束条件:
∑
i
=
1
m
a
i
?
l
a
b
e
l
(
i
)
=
0
\sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
i=1∑m?ai??label(i)=0
3、完整示例代码
from operator import mul
import numpy as np
from numpy import *
from numpy.core.defchararray import multiply
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
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
return j
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose()
b = 0
m, n = shape(dataMatrix)
alphas = np.mat(zeros((m, 1)))
iter = 0
while iter < maxIter:
alphaPairsChange = 0
for i in range(m):
fXi = float(np.multiply(alphas, labelMat).T *
(dataMatrix*dataMatrix[i, :].T))+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 *
(dataMatrix*dataMatrix[j, :].T))+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, C+alphas[j]+alphas[i])
if L == H:
print("L==H")
continue
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
alphas[j] -= labelMat[j]*(Ei-Ej)/eta
alphas[j] = clipAlpha(alphas[j], H, L)
if abs(alphas[j]-alphaJold) < 0.00001:
print("j not moving enough")
continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
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
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
alphaPairsChange += 1
print("Iter: %d i: %d, pairs changed %d" %
(iter, i, alphaPairsChange))
if alphaPairsChange == 0:
iter += 1
else:
iter = 0
print("iteration number: %d" % iter)
return b, alphas
dataArr, labelArr = loadDataSet('testSet.txt')
smoSimple(dataArr, labelArr, 0.6, 0.001, 10)
利用完整Platt SMO算法加速优化
??在几百个点组成的小规模数据集上,简化版SMO算法的运行是没有问题的,但是在更大的数据集上运行速度就会变慢。完整的Platt SMO算法唯一不同的就是选择alpha的方式。完整版Platt SMO算法应用了一些能够提高启发的方法。 ??Platt SMO算法是通过一个外循环来选择第一个alpha值的,并且其选择过程会在两种之间进行交替:一种是在所有数据集上进行单遍扫描,另一种方式是在非边界alpha中实现单遍扫描,所谓非边界alpha指的是那些不等于边界的0或C的alpha值。对整个数据集的扫描比较容易,而实现非边界alpha值的扫描,首先需要建立这些alpha值的列表,然后在对这个表进行遍历。同时该步骤还会跳过那些已知不会改变的alpha的值。 ??在选择第一个alpha值后,算法会通过一个内循环来选择第二个alpha值,在优化过程中会通过最大化步长的方式来获取第二个alpha值。
完整示例代码
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
class optStruct:
def __init__(self, dataSetIn, classLabels, C, toler):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
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
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 selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.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)
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]
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
return j
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
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)
alphaIoid = oS.alphas[i].copy()
alphaJoid = 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.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 0
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
updateEk(oS, i)
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.X[i, :]*oS.X[i,
:].T-oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.X[i, :]*oS.X[j, :].T
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.X[i, :]*oS.X[j,
:].T-oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*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:
b = (b1+b2)/2.0
return 1
else:
return 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 changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
else:
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, pair changed %d" %
(iter, i, alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
print("Iteration number: %d" % iter)
return oS.b, oS.alphas
def calcWs(alphas, dataArr, classLabels):
X = np.mat(dataArr)
labelMat = np.mat(classLabels).transpose()
m, n = shape(X)
w = zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i], X[i, :].T)
return w
def plotFit(dataMat, labelMat, alphas, b):
dataArr = np.array(dataMat)
n = np.shape(dataArr)[0]
xcord1, ycord1 = [], []
xcord2, ycord2 = [], []
for i in range(n):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i, 0])
ycord1.append(dataArr[i, 1])
else:
xcord2.append(dataArr[i, 0])
ycord2.append(dataArr[i, 1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-1.0, 8.0, 0.1)
weights = np.multiply(alphas, np.mat(
labelMat).transpose()).T * np.mat(dataMat)
weights = np.array(weights)[0]
y = (-np.array(b)[0][0]-weights[0]*x)/weights[1]
ax.plot(x, y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
dataArr, labelArr = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
plotFit(dataArr, labelArr, alphas, b)
完
整
版
P
l
a
t
t
S
M
O
算
法
结
果
完整版Platt SMO算法结果
完整版PlattSMO算法结果
在复杂数据上应用核函数
??核函数(kernel)目的是:将数据转换成易于分类器理解的形式。
1、利用核函数将数据映射到高维空间
??从一个特征空间到另一个特征空间的映射,这种映射会将低维特征空间映射到高维空间。这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。核函数相当于一个包装器或者接口,它能把数据从某个很难处理的形式转换成一个较容易处理的形式。 ??SVM优化中一个特别好的地方就是所有的运算都可以写成内积(也称点积)的形式。向量的内积是两个向量相乘,之后得到单个标量或者数值。把内积运算换成核函数的方式称为核技巧或者核变化
2、径向基核函数
??径向基核函数是SVM中常用的一个核函数。径向基核函数是一个采用向量作为自变量的函数,能够基于向量距离运算输出一个标量。这个距离可以是从<0,0>向量或者其他向量开始计算的距离。径向基核函数的高斯版本:
k
(
x
,
y
)
=
e
x
p
?
∥
x
?
y
∥
2
2
σ
2
?
?
k(x,y)=exp\lgroup\dfrac {{\parallel x-y\parallel}^2}{2\sigma^2} \rgroup?
k(x,y)=exp?2σ2∥x?y∥2??? 其中,
σ
\sigma
σ是用户定义的确定达到达率或者说函数值跌落到0的速度参数。 ??上述高斯核函数将数据从其特征空间映射到更高维的空间,具体说是映射到一个无穷维的空间。
3、测试中使用核函数
??构建一个对数据点进行有效分类的分类器该分类器使用了径向基核函数,径向基核函数有一个用户定义的输入
σ
\sigma
σ,我么需要先确定它的大小,再利用该核函数构建一个分类器。
4、完整示例代码
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
def kernelTrans(X, A, kTup):
m, n = shape(X)
K = np.mat(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))
else:
raise NameError(
'Houston We Have a Problem -- That Kernel is not recognized')
return K
class optStruct:
def __init__(self, dataSetIn, classLabels, C, toler, kTup):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
self.K = np.mat(zeros((self.m, self.m)))
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
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)
alphaIoid = oS.alphas[i].copy()
alphaJoid = 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.C+oS.alphas[j]+oS.alphas[i])
if L == H:
print("L==H")
return 0
eta = 2.0*oS.K[i, j]-oS.K[i, i]-oS.K[j, j]
if eta >= 0:
print("eta>=0")
return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
updateEk(oS, i)
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, i] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[i, j]
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, j] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*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:
b = (b1+b2)/2.0
return 1
else:
return 0
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
def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.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)
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]
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
return j
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
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
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
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):
alphaPairsChanged = 0
if entireSet:
for i in range(oS.m):
alphaPairsChanged += innerL(i, oS)
print("FullSet, iter: %d i: %d, pairs changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
else:
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, pair changed %d" %
(iter, i, alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
print("Iteration number: %d" % iter)
return oS.b, oS.alphas
def calcWs(alphas, dataArr, classLabels):
X = np.mat(dataArr)
labelMat = np.mat(classLabels).transpose()
m, n = shape(X)
w = zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i], X[i, :].T)
return w
def testRbf(k1=1.3):
dataArr, labelArr = loadDataSet('testSetRBF.txt')
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
dataMat = np.mat(dataArr)
labelMat = np.mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]
sVs = dataMat[svInd]
labelSV = labelMat[svInd]
print("There are %d Support Vectors" % shape(sVs)[0])
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], ('rbf', k1))
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
print("The Train error rate is: %f" % (float(errorCount)/m))
dataArr, labelArr = loadDataSet('testSetRBF2.txt')
errorCount = 0
dataMat = np.mat(dataArr)
labelMat = np.mat(labelArr).transpose()
m, n = shape(dataMat)
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], ('rbf', k1))
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
print("The test error rate is: %f" % (float(errorCount)/m))
testRbf()
示例:手写识别问题回顾
基于SVM的数字识别一般流程
- 收集数据:提供的文本文件
- 准备数据:基于二值图像的构造向量
- 分析数据:对图像向量进行目测
- 训练算法:采用两种不同的核函数,并对径向基核函数采用不同的设置来运行SMO算法
- 测试算法:编写一个函数来测试不同的核函数并计算错误率
- 使用算法:一个图像识别的完整应用还需要一些图像处理的知识
??第二章时,介绍了KNN算法来识别数字,但是需要保留样本,采用SVM其需要的样本数量少了很多。
完整示例代码展示
from os import listdir
import numpy as np
from numpy import *
def loadImages(dirName):
hwLabels = []
trainingFileList = listdir(dirName)
m = len(trainingFileList)
trainingMat = 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
class optStruct:
def __init__(self, dataSetIn, classLabels, C, toler, kTup):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
self.K = np.mat(zeros((self.m, self.m)))
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
def kernelTrans(X, A, kTup):
m, n = shape(X)
K = np.mat(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))
else:
raise NameError(
'Houston We Have a Problem -- That Kernel is not recognized')
return K
def img2vector(filename):
resultVect = np.zeros((1, 1024))
fr = open(filename)
for i in range(32):
linesStr = fr.readline()
for j in range(32):
resultVect[0, 32*i+j] = int(linesStr[j])
return resultVect
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
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):
alphaPairsChanged = 0
if entireSet:
for i in range(oS.m):
alphaPairsChanged += innerL(i, oS)
print("FullSet, iter: %d i: %d, pairs changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
else:
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, pair changed %d" %
(iter, i, alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
print("Iteration number: %d" % iter)
return oS.b, oS.alphas
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)
alphaIoid = oS.alphas[i].copy()
alphaJoid = 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.C+oS.alphas[j]+oS.alphas[i])
if L == H:
print("L==H")
return 0
eta = 2.0*oS.K[i, j]-oS.K[i, i]-oS.K[j, j]
if eta >= 0:
print("eta>=0")
return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
updateEk(oS, i)
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, i] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[i, j]
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, j] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*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:
b = (b1+b2)/2.0
return 1
else:
return 0
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
def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.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)
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]
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
return j
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadImages('../机器学习02/trainingDigits')
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
dataMat = np.mat(dataArr)
labelMat = np.mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]
sVs = dataMat[svInd]
labelSV = labelMat[svInd]
print("There are %d Support Vectors" % shape(sVs)[0])
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
print("The Train error rate is: %f" % (float(errorCount)/m))
dataArr, labelArr = loadImages('../机器学习02/testDigits')
errorCount = 0
dataMat = np.mat(dataArr)
labelMat = np.mat(labelArr).transpose()
m, n = shape(dataMat)
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
print("The test error rate is: %f" % (float(errorCount)/m))
testDigits(('rbf', 20))
提示:这个代码和上面的代码几乎一样,只是个别地方做了修改
小结
??支持向量机是一种分类器。之所以称为 “机” 是因为它会产生一个二值决策结果,即它是一种决策 “机” 。支持向量机的泛化错误率较低,具有良好的学习能力,且学到的结果具有很好的推广性。这些优点使支持向量机十分流行。 ??支持向量机通过求解一个二次优化问题来最大化分类化间隔。引入SMO算法加快了SVM的训练速度。核方法或者是核技巧将数据(优势是非线性的数据)从一个低维空间映射到高维空间,可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题来求解。 ??支持向量机是一个二类分类器。当用其解决多类问题时,需要额外的方法对其进行扩展。SVM的效果也对优化参数和所用核函数中的参数敏感。
提醒: 本节的数据较多,不再展示,有需要的话,请自行查找。
|