11 min to read
支持向量机原理及利用python实现
- content {:toc}
朴素贝叶斯原理
原理
定义超平面$h(x)=\omega^Tx+b$，样本$x$到超平面的几何间隔$r=\frac{\omega^Tx+b}{||\omega||}$，其中 $||\omega||=\sqrt{\omega_0^2+\omega_1^2+...+\omega_n^2}$。
在支持向量机中把几何间隔$r$作为最大化间隔进行分析，并采用-1和1作为类别标签。我们把其中一个支撑向量$x^*$到最优超平面的距离定义为：
$$ r^=\frac{h(x^)}{||\omega||}= \begin{cases} \frac{1}{||\omega||}& if: y^=h(x^)=+1\ -\frac{1}{||\omega||}& if: y^=h(x^)=-1 \end{cases} $$
对于任意空间点$x_i$，有约束条件：
$$ \begin{cases} \omega^Tx_i+b\geq1& y_i=+1\ \omega^Tx_i+b\leq-1& y_i=-1 \end{cases} $$
优化问题变为：
$$ \begin{cases} \mathop{max}\limits_{\omega,b}\frac{2}{||\omega||} \iff \mathop{min}\limits_{\omega,b}\frac{1}{2}||\omega||^2\ y_i(\omega^Tx_i+b)\geq1,(i=1,2,...,n) \end{cases} $$
拉格朗日乘子发对偶问题
上式的优化目标是二次的，约束是线性的，整体是一个凸二次规划问题。将其转化为拉格朗日对偶问题后更容易求解，也方便后面引入核函数。
$$ L(x)=f(x)+\sum\alpha g(x) $$
其中，$f(x)$是目标函数，$g(x)$是约束条件。为了得到拉格朗日函数的最优解。
有n个点设置对应的拉格朗日乘子$\alpha_i\ge 0$，得到原始问题的拉格朗日函数：
$$ L(\omega,b,\alpha)=\frac{1}{2}||\omega||^2+\sum_{i=1}^n\alpha_i[1-y_i(\omega^Tx_i+b)] $$
目标是让$L(\omega,b,\alpha)$针对$\alpha$达到最大值，然后再让$L(\omega,b,\alpha)$针对$\omega,b$达到最小值。目标函数
$$ \mathop{min}\limits_{\omega,b}\mathop{max}\limits_{\alpha}L(\omega,b,\alpha) $$
其对偶问题是：
$$ \mathop{max}\limits_{\alpha}\mathop{min}\limits_{\omega,b}L(\omega,b,\alpha) $$
需满足下列KKT条件才能有解：
$$ \begin{cases} \alpha_i\geq0\ y_i(\omega^Tx_i+b)-1\geq0\ \alpha_i[y_i(\omega^Tx+b)-1]=0 \end{cases} $$
$\omega,b$的部分
$L(\omega,b,\alpha)$对$\omega$和$b$分别求偏导，并令其等于0，带入原拉格朗日式子得到对偶问题：
$$ \begin{cases} \mathop{max}\limits_{\alpha}\sum_{i=1}^{n}\alpha_i-\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_jx_i^Tx_j\ s.t.\qquad\sum_{i=1}^{n}\alpha_iy_i=0\ \alpha_i\geq0\ \end{cases} $$
$\alpha$的部分
利用SMO算法求解，就是将最大优化问题分解为多个小优化问题来求解，这些小优化问题往往很容易求解，并且对他们进行顺序求解的结果与将他们作为一个整体求解的结果完全一致。
软间隔
引入松弛变量，允许模型出现某些错误。
$$ \begin{cases} \mathop{min}\limits_{\omega,b}\frac{1}{2}|\omega||^2+C\sum_{i=1}^n\xi_i\ y_i(\omega^Tx_i+b)\geq1-\xi_i,\xi_i\geq0,(i=1,2,...,n) \end{cases} $$
得到新的拉格朗日函数
$$ L(\omega,b,\xi,\alpha,\mu)=\frac{1}{2}||\omega||^2+C\sum_{i=1}^n\xi_i+\sum_{i=1}^n\alpha_i[1-\xi_i-y_i(\omega^Tx_i+b)]-\sum_{i=1}^n\mu_i\xi_i $$
需满足下列KKT条件才能有解：
$$ \begin{cases} \alpha_i\geq0\,\mu_i\geq0 y_i(\omega^Tx_1+b)-1+\xi_i\geq0\ \alpha_i[y_i(\omega^Tx_i+b)-1+\xi_i]=0\ \xi_i\geq0,\mu_i\xi_i=0 \end{cases} $$
$L(\omega,b,\xi,\alpha,\mu)$对$\omega,b和\xi$分别求偏导，并令其等于0，带入原拉格朗日式子得到对偶问题：
$$ \begin{cases} \mathop{max}\limits_{\alpha}\sum_{i=1}^{n}\alpha_i-\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_jx_i^Tx_j\ s.t.\qquad\sum_{i=1}^{n}\alpha_iy_i=0\ 0\leq\alpha_i\leq C \end{cases} $$
核函数
在现实任务中，可能并不存在一个超平面将数据集完美地线性分割开，这样需要将原始空间映射到一个更高维的空间，如果在这一高伟空间数据集变得线性可分，那么问题就得到了解决。
超平面$\omega^T{\color{red} {\phi(x)}}+b=0$
经过推导得到：
$$ \begin{cases} \mathop{max}\limits_{\alpha}\sum_{i=1}^{n}\alpha_i-\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j{\color{red} {\phi(x_i)^T\phi(x_j)}}\ s.t.\qquad \sum_{i=1}^{n}\alpha_iy_i=0\ \alpha_i\geq0\ \end{cases} $$
如果存在一个函数$\kappa(x_i,x_j)=\phi(x_i)^T\phi(x_j)$就可以在低维空间计算高维空间的点积结果。称$\kappa(x_i,x_j)$为核函数。
常见的核函数有：
| 名称 | 表达式 | 备注 |
|---|---|---|
| 线性核 | $\kappa(x_i,x_j)=x_i^Tx_j$ | |
| 多项式核 | $\kappa(x_i,x_j)=(x_i^Tx_j)^d$ | $d\geq1$为多项式次数 |
| 高斯核 | $\kappa(x_i,x_j)=exp(-\frac{\vert\vert x_i-x_j\vert\vert^2}{2\sigma^2}$ | $\sigma>0$为高斯核宽带 |
| 拉普拉斯核 | $\kappa(x_i,x_j)=exp(-\frac{\vert\vert x_i-x_j\vert\vert}{\sigma}$ | $\sigma>0$ |
| Sigmoid核 | $\kappa(x_i,x_j)=tanh(\beta x_i^Tx_j+\theta)$ | $\beta>0,\theta<0,tanh$为双曲正切函数 |
实战
所需数据
代码
import matplotlib.pyplot as plt
import numpy as np
import random
"""
函数说明:读取数据
Parameters:
fileName - 文件名
Returns:
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
"""
函数说明:随机选择alpha
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def selectJrand(i, m):
j = i #选择一个不等于i的j
while (j == i):
j = int(random.uniform(0, m))
return j
"""
函数说明:修剪alpha
Parameters:
aj - alpha_j值
H - alpha上限
L - alpha下限
Returns:
aj - alpah值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
函数说明:数据可视化
Parameters:
dataMat - 数据矩阵
labelMat - 数据标签
Returns:
无
"""
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()
"""
函数说明:简化版SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
无
"""
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
"""
函数说明:分类结果可视化
Parameters:
dataMat - 数据矩阵
w - 直线法向量
b - 直线解决
Returns:
无
"""
def showClassifer(dataMat, w, b):
#绘制样本点
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], s=30, alpha=0.7) #正样本散点图
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图
#绘制直线
x1 = max(dataMat)[0]
x2 = min(dataMat)[0]
a1, a2 = w
b = float(b)
a1 = float(a1[0])
a2 = float(a2[0])
y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
plt.plot([x1, x2], [y1, y2])
#找出支持向量点
for i, alpha in enumerate(alphas):
if abs(alpha) > 0:
x, y = dataMat[i]
plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
plt.show()
"""
函数说明:计算w
Parameters:
dataMat - 数据矩阵
labelMat - 数据标签
alphas - alphas值
Returns:
无
"""
def get_w(dataMat, labelMat, alphas):
alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
return w.tolist()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet('testSet.txt')
b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
w = get_w(dataMat, labelMat, alphas)
showClassifer(dataMat, w, b)
参考资料
[1]. SVM原理以及Tensorflow 实现SVM分类(附代码)
[2]. 【机器学习】支持向量机SVM原理及推导
Comments