SVM编程实现python

深入解析python版SVM源码系列--简化版SMO算法

SVM使用SMO算法来解决其中涉及到的二次规划问题。一个简单版本的SMO算法的实现如下：

''' 随机选择随机数，不等于J '''
def selectJrand(i,m):
j=i #we want to select any J not equal to i
while (j==i):
j = int(random.uniform(0,m))  # 一直在挑选随机数j，直到不等于i，随机数的范围在0~m
return j  # 返回挑选好的随机数

''' 门限函数 '''
def clipAlpha(aj,H,L):  # 最大不能超过H，最小不能低于L
if aj > H:
aj = H
if L > aj:
aj = L
return aj

''' 简化版的SMO函数 '''
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):  # 输入数据，标记，常数C，容错率，最大迭代次数
dataMatrix = mat(dataMatIn);   # 转换成矩阵
labelMat = mat(classLabels).transpose()  # 转换成矩阵，并转置，标记成为一个列向量，每一行和数据矩阵对应
m,n = shape(dataMatrix)  # 行，列

b = 0;  # 参数b的初始化
alphas = mat(zeros((m,1)))  # 参数alphas是个list，初始化也是全0，大小等于样本数
iter = 0  # 当前迭代次数，maxIter是最大迭代次数

while (iter < maxIter):  # 当超过最大迭代次数，推出
alphaPairsChanged = 0  # 标记位，记录alpha在该次循环中，有没有优化
for i in range(m):  # 第i个样本
fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b  # 第i样本的预测类别
Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions # 误差

#是否可以继续优化
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m)  # 随机选择第j个样本
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b  # 样本j的预测类别
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 * 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)  # 门限函数阻止alpha_j的修改量过大

#如果修改量很微小
if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue

# alpha_i的修改方向相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
#the update is in the oppostie direction
# 为两个alpha设置常数项b
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

# 说明alpha已经发生改变
alphaPairsChanged += 1
print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)

#如果没有更新，那么继续迭代；如果有更新，那么迭代次数归0，继续优化
if (alphaPairsChanged == 0): iter += 1
else: iter = 0
print "iteration number: %d" % iter

# 只有当某次优化更新达到了最大迭代次数，这个时候才返回优化之后的alpha和b
return b,alphas

使用SMO算法可以迭代更新参数alphas,而依据alphas我们可以得到最终的分类决策超平面的参数w和b。

对于SMO算法中有一个重要的代码：计算样本的预测类别。如下：

fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b  # 第i样本的预测类别

我们知道原始的预测类别计算公式是用决策面的参数w和b表示的，那么为什么这里的貌似不一样呢？

原始的预测类别计算公式为：



其中w可以表示为：



然后分类函数可以转化为：



关于这个的解释，july博客上说的比较清晰：

这里的形式的有趣之处在于，对于新点 x的预测，只需要计算它与训练数据点的内积即可（<.>表示向量内积），这一点至关重要，是之后使用 Kernel 进行非线性推广的基本前提。

这样子的表示形式和上面代码就一致了。

这儿还有一个现象可以分析出来：哪些是支持向量。

答：alpha不等于0的为支持向量。

所谓 Supporting Vector 也在这里显示出来??事实上，所有非Supporting Vector 所对应的系数alpha都是等于零的，因此对于新点的内积计算实际上只要针对少量的“支持向量”而不是所有的训练数据即可。

为什么非支持向量对应的alpha等于零呢？直观上来理解的话，就是这些“后方”的点??正如我们之前分析过的一样，对超平面是没有影响的，由于分类完全有超平面决定，所以这些无关的点并不会参与分类问题的计算，因而也就不会产生任何影响了。



注意到如果 xi 是支持向量的话，上式中红颜色的部分是等于 0 的（因为支持向量的 functional margin 等于 1 ），而对于非支持向量来说，functional margin 会大于 1 ，因此红颜色部分是大于零的，而alpha_i又是非负的，为了满足最大化，必须alpha_i等于 0 。这也就是这些非Supporting Vector 的点的局限性。

所以，在我们运行SMO算法程序之后，可以根据这个特点求得支持向量，也就是alpha不等于0。



SMO算法在SVM源码中的更新步骤是最为重要的，如下：

#是否可以继续优化
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m)  # 随机选择第j个样本
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b  # 样本j的预测类别
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 * 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)  # 门限函数阻止alpha_j的修改量过大

#如果修改量很微小
if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue

# alpha_i的修改方向相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
#the update is in the oppostie direction
# 为两个alpha设置常数项b
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

# 说明alpha已经发生改变
alphaPairsChanged += 1
print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)

这个算法步骤是这样子的：

Step 1：首先随机选择两个alpha_i和alpha_j，如果它们满足：



说明它们可以进行优化，也就是不满足KKT条件。

对应的源码部分是：

if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):

Step 2：如果选择的两个alpha是可以更新优化的，那么我们使用下面的算法（SMO的核心）来学习两个新的alpha：