SVM 的思想

有时候,看似不可分的数据,我们换一个视角或者维度,它就可分了:

ejSVM

为了可分,必须增加额外维度,然后找到一个最优分离超平面

对于前者,要采用不同的核函数。对于后者,如何衡量最优?引入了间隔的概念。

线性可分的情况

即能够存在平面,使得两侧恰为不同类的样本。

间隔 / 边距:即下面的灰色区域。定义间隔为 $\rho$.

最大间隔超平面:即下面灰色区域正中间的黑线。

支持向量:即下面落在灰色边界上的向量。

image-20211129204439476

为了找到最合适的分界线,就相当于要找到最大间隔超平面。

超平面的表达式为 $X^TW +b$ (或者 $\mathbf {w} \cdot \mathbf {x} +b=0$,这要求 $X,T$ 都是一维向量):

  • $X$:参数向量,可以变
  • $W$:法向量。规定法向量的正向为正例类。也称权重向量。
  • $b$:截距。

如何表示点 $x_i$ 距离超平面的远近?

根据点到直线的距离公式,距离为

$$ d = \frac{1}{\|w\|}|\mathbf w \cdot \mathbf {x_0} + b| $$

其中 $||w||$​ 是欧几里得范数(其实就是向量间距离)。

范数:详见本人 数学分析笔记 / R 拓扑学基础

什么是函数间隔和几何间隔?

一般来说,一个点距离分离超平面的远近,反应分类预测确信度的高低。在超平面 $\mathbf {w} \cdot \mathbf {x}+b=0$ 确定的情况下,$|\mathbf {w} \cdot \mathbf {x} +b|$ 能够相对地表示向量点 $x_0$ 距离超平面的远近。而 $\mathbf {w} \cdot \mathbf {x_0} +b$ 的符号与类标签 $y$ 的符号($\pm1$)是否一致能够表示分类是否正确。所以可用量 $y (\mathbf {w} \cdot \mathbf {x} +b)$ 来表示分类的正确性及确信度,这就是函数间隔(functional margin)的概念。

超平面关于样本点的函数间隔

对于给定的训练数据集 $T$ 和超平面 $(w,b)$,定义超平面 $(w,b)$ 关于样本点 $(x_i,y_i)$ 的函数间隔为

$$ \hat{\gamma}_i=y_i(\mathbf{w} \cdot \mathbf{x_i} +b) $$

超平面对于训练集的函数间隔

定义超平面关于训练数据集 $T$ 的函数间隔为超平面 $(w,b)$ 关于 T 中所有样本点 $(x_i,y_i)$ 的函数间隔之最小值,即

$$ \hat{\gamma}=\mathop{\text{min}}_{i=1,...,N}\hat{\gamma}_i $$

image-20211129204439476

为什么定义为最小值?上图体会

规范化 2 - 范数,得到几何间隔

只要成比例地改变 $w$ 和 $b$,例如将它们改为 $2w$ 和 $2b$,超平面并没有改变,但函数间隔却成为原来的 $2$ 倍,这样不但造成算力浪费,在数学上也是不优雅的。因此,我们让 $\mathbf {w}$ 和 $b$ 都乘以 $\mathbf {w}$ 的长度,从而固定下来:

$$ {\gamma}_{i}=y_{i}(\frac{\mathbf{w} }{||w||} \mathbf{x_{i}}+\frac{b}{||w||}) $$

这称为几何间隔。用以表示样本点 $(\mathbf {x_i}, y_i)$ 到超平面的距离。

上式说明了几何间隔和函数间隔的关系:$\gamma _i = \dfrac {\hat {\gamma _i}}{||w||}$ 当 $||w|| = 1$ 时,几何间隔就是函数间隔

image-20211129204439476

再看一次,我们要得到正确的间隔(灰色区域),就得取几何间隔的最小值。

因此:定义超平面 $(w,b)$ 关于训练数据集 $T$ 的几何间隔为超平面 $(w,b)$ 关于 $T$ 所有样本点 $(x_i,y_i)$ 的几何间隔之最小值,即

$$ {\gamma}=\min_{i=1,...,N}{\gamma}_{i} $$

超平面 $(w,b)$ 关于样本点 $(x_i,y_i)$ 的几何间隔是实例点到超平面的带符号距离(signed distance),当样本点被超平面正确分类时就是实例点到超平面的距离。

为啥带符号?因为要区分正例和反例,所以利用在平面两侧区分符号不同这一特点区分

好,那么问题来了,我们想找到最好的平面,而最好的平面,要满足既是最小、也是最大。

  • 最小:取各个样本点的几何距离,其中的最小值作为 $\gamma $
  • 最大:要在满足上面这个条件的情况下,找到最大的 $\gamma $

也即目标函数:

$$ \max \gamma = \max \min \gamma _i $$

得到满足时的平面 $(\mathbf {w}, b)$ 。

怎么找?

如何找最大间隔时平面参数?

这步操作称为 硬间隔最大化(线性可分,所以是间隔)。

找最大硬间隔,就是找最优化问题:

$$ \begin{aligned} &\mathop{\text{max}}_{\mathbf{w},b}\quad \gamma\\ & \text{ s.t.} \quad \gamma _i \geqslant \gamma,\ i=1,2,...,N\\ \end{aligned} $$

的解。

$\mathop {\text {max}}_{\mathbf {w},b} \gamma $ 表示我们想要使 $\gamma $ 最大化,然后得到此时的 $\mathbf {w}, b$ 其中的 $\gamma _i = y_{i}\left (\frac { \mathbf {w} }{\|w\|} x_{i}+\frac {b}{\|w\|}\right)$

代入几何间隔和函数间隔的关系($\gamma _i = \dfrac {\hat {\gamma _i}}{||w||}$ ),目标函数变形为:

$$ \mathop{\text{max}}_{\mathbf{w},b}\quad \dfrac{\hat{\gamma }}{||w||} $$

另外再代入一次这个关系到 $\text {s.t.}$ 那儿,得到:

$$ \text{ s.t. } \quad y_i(\mathbf{w} \mathbf{x_i} + b) \geqslant \hat{\gamma } $$

每一个初中生都知道,如果 $|x|$ 变大,那么 $\dfrac {1}{|x|}$ 变小。所以我们可以把目标函数改写成:

$$ \mathop{\text{min}}_{\mathbf{w},b}\quad \dfrac{||w||}{\hat{\gamma }} $$

同时我们也说过,就如你把 $x + y - 1 = 0$ 写成 $2x + 2y -2 = 0$ ,函数间隔可以随便取。取一个 1. 目标函数变成:

$$ \mathop{\text{min}}_{\mathbf{w},b}\quad ||w|| $$

为了好算,反正只要不影响最后 $\mathbf {w}, b$ 的正确性,把目标函数再次变成:

$$ \mathop{\text{min}}_{\mathbf{w},b}\quad \frac{1}{2}||w||^2 $$

为什么可以加个平方?因为 $y = x ^{2}$ 在 $[0, + \infty)$ 单增。 为什么不是 1 次方或 3 次方?因为有人试过,不好算。 为什么会有个 $1/2$ ?因为后面要求导,方便。

$$ \begin{aligned} &\mathop{\text{min}}_{\mathbf{w},b}\quad \frac{1}{2}||w||^2\\ &\text{s.t.} \quad {\color{#228a96}{ y_i\left(\mathbf{w} \cdot \mathbf{x_i} +b\right)-1\geqslant 0,\ i=1,2,\cdots ,N}} \\ \end{aligned} $$

为什么会有蓝色这行?你搜前文 “再代入一次这个关系” 这句话,然后代入函数间隔 $\hat {\gamma} = 1$ 即可得到。

如果求出了约束最优化问题的解 $\mathbf {w}, b$,那么就可以得到最大间隔分离超平面及分类决策函数 $f (x)=\operatorname {sign} (\mathbf {w} \mathbf {x} + b)$,即线性可分支持向量机模型

问题来了,怎么解呢?

如何求解最优化问题?

本文假定你已经了解过 拉格朗日乘子法。即通过极值处的梯度共线条件($\nabla f = \lambda \nabla g$)列方程组求解极值。

对于 SVM 的最优化问题,限制条件( s.t. 后面那个)实际上有一堆($i = 1, 2, \cdots ,n$ )所以我们引入拉格朗日算子 $\alpha _i \mid_{i = 1, 2, \cdots , n}$,$\alpha=(\alpha_1,\alpha_2,…,\alpha_n)$ 为拉格朗日乘子向量。

设 Lagrange 函数为:

$$ L(w,b,\alpha)=\frac{1}{2}||w||^2-\sum_{i=1}^n\alpha_iy_i(w\cdot x_i+b)+\sum_{i=1}^n\alpha_i $$

则原问题转化为:

$$ \begin{aligned} &\mathop{\text{min}}_{w,b}\mathop{\text{max}}_{\alpha}\quad L(w,b,\alpha)\\ &\text{s.t.}\quad \alpha_i\geqslant 0\\ \end{aligned} $$

注意,其中消去了 条件 $1-y_i (w^Tx_i+b)\leqslant 0$,之所以可以这么做: 因为 $ y*i\left (\mathbf {w} \cdot \mathbf {x_i} +b\right)-1\geqslant 0$ 时, 成立 $\mathop {\text {max}}*{\alpha} L (w,b,\alpha)=\frac {1}{2}||w||^2+0=\frac {1}{2}||w||^2$ 而 $ y_i\left (\mathbf {w} \cdot \mathbf {x_i} +b\right)-1< 0$ 时,$\mathop {\text {max}}_{\alpha} L (w,b,\alpha)=\frac {1}{2}||w||^2+\infty=\infty$ 即不存在解。

根据拉格朗日对偶性,原始问题的对偶问题是:

可以参考 拉格朗日对偶(Lagrange duality)、KKT 条件,但我们暂时大可不必深究,否则这门课学不完的。参考书 Convex Optimization – Boyd and Vandenberghe.

$$ \mathop{\text{max}}_{\alpha}\mathop{\text{min}}_{w,b}L(w,b,\alpha) $$

这里的交换严格来说需要证明前后等价,这涉及到强对偶、若对偶、Slater 条件云云,暂时别管那么多。毕竟咱们只是本科生,一节课还学这么多东西,凸优化也没学过呜呜。

求解 $\mathop {\text {min}}_{w,b} L (w,b,\alpha)$

$L (w,b,\alpha)$ 分别对 $w,b$ 求偏导数,并令其等于 0。

$$ \begin{aligned} \bigtriangledown_wL(w,b,\alpha)&=w-\sum_{i=1}^N\alpha_iy_ix_i=0\\ \bigtriangledown_bL(w,b,\alpha)&=-\sum_{i=1}^N\alpha_iy_i=0\\ \end{aligned} $$

$$ \begin{aligned} &w=\sum_{i=1}^N\alpha_iy_ix_i\\ &\sum_{i=1}^N\alpha_iy_i=0\\ \end{aligned} $$

从而

$$ \begin{aligned} L(w,b,\alpha)&=\frac{1}{2}||w||^2-\sum_{i=1}^N\alpha_iy_i(w\cdot x_i+b)+\sum_{i=1}^N\alpha_i\\ &=\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(\mathbf{x_i}\cdot \mathbf{x_j})-\sum_{i=1}^N\alpha_iy_i\left(\left( \sum_{j=1}^N\alpha_jy_jx_j \right)\cdot x_i+b\right)+\sum_{i=1}^N\alpha_i\\ &=-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(\mathbf{x_i}\cdot \mathbf{x_j})+\sum_{i=1}^N\alpha_i\\ \end{aligned} $$

即:

$$ \mathop{\text{min}}_{w,b}L(w,b,\alpha)=-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(\mathbf{x_i}\cdot \mathbf{x_j})+\sum_{i=1}^N\alpha_i $$

现在问题变成:

$$ \begin{aligned} &\mathop{\text{max}}_{\alpha}\ -\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(\mathbf{x_i}\cdot \mathbf{x_j})+\sum_{i=1}^N\alpha_i\\ & \text{ s.t. } \quad \sum_{i=1}^N\alpha_iy_i=0\\ &\quad\quad\quad\alpha_i\geqslant 0,\quad i=1,2,...,N \end{aligned} $$

我们取个相反数,将最大值问题转换成最小值问题:

$$ \begin{aligned} &\mathop{\text{min}}_{\alpha}\ \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(\mathbf{x_i}\cdot \mathbf{x_j})-\sum_{i=1}^N\alpha_i\\ & \text{ s.t. } \quad \sum_{i=1}^N\alpha_iy_i=0\\ &\quad\quad\quad\alpha_i\geqslant 0,\quad i=1,2,...,N \end{aligned} $$

上述目标函数假设样本点线性硬可分。

同时解决线性不可分

但现实世界并非完美。线性不可分意味着某些样本点不能满足线性可分中的函数间隔大于等于 1 的约束条件。

为了提高容错性,为每个样本点 $(\mathbf {x}_i, y_i)$ 引入 松弛变量 $ξ_i$,使函数间隔加上松弛变量大于等于 1。这样,约束条件变为:

$$ y_i(\mathbf{w} \cdot \mathbf{x} _i+b) \geqslant 1-\xi_i $$

为啥要这么引入?等我学完凸优化再来回答吧……

再乘以惩罚因子 $C$,从而目标函数变为

$$ \mathop{\text{min}}_{\mathbf{w},b}\quad \frac{1}{2}\|w\|^{2}+C \sum_{i=1}^{N} \xi_{i} $$

优化问题变成:

$$ \begin{aligned} &\mathop{\text{min}}_{\mathbf{w},b,\xi}\quad \frac{1}{2}\|w\|^{2}+C \sum_{i=1}^{N} \xi_{i}\\ &\text{s.t.}\ \ \quad y_i(w\cdot x_i+b)\geqslant 1-\xi_i,\ i=1,2,...,N\\ &\ \ \ \quad \quad \xi_i\geqslant 0,\ i=1,2,...,N\\ \end{aligned} $$

最优化问题的 Lagrange 函数变为:

$$ \begin{align} &L(w, b, \xi, \alpha, \mu)=\\ &\quad \frac{1}{2}\|w\|^{2}+C \sum_{i=1}^{N} \xi_{i}-\sum_{i=1}^{N} \alpha_{i}\left(y_{i}\left(w \cdot x_{i}+b\right)-1+\xi_{i}\right)-\sum_{i=1}^{N} \mu_{i} \xi_{i} \end{align} $$

通过拉格朗日函数我们可以将原始问题

$$ \begin{aligned} &\mathop{\text{min}}_{w,b,\xi}\quad \frac{1}{2}||w||^2+C\sum_{i=1}^N\xi_i\\ &\text{s.t.}\ \ \quad y_i(w\cdot x_i+b)\geqslant 1-\xi_i,\ i=1,2,...,N\\ &\ \ \ \quad \quad \xi_i\geqslant 0,\ i=1,2,...,N\\ \end{aligned} $$

转化为:

$$ \begin{aligned} &\mathop{\text{min}}_{w,b,\xi}\mathop{\text{max}}_{\alpha,\mu}\quad L(w,b,\xi,\alpha,\mu)\\ &\text{s.t.}\quad \alpha_i\geqslant 0, \mu_i\geqslant0\\ \end{aligned} $$

对偶问题:

$$ \mathop{\text{max}}_{\alpha,\mu}\mathop{\text{min}}_{w,b,\xi}\quad L(w,b,\xi,\alpha,\mu) $$

$$ \begin{aligned} &\bigtriangledown_wL(w,b,\xi,\alpha,\mu)=w-\sum_{i=1}^N\alpha_iy_ix_i=0\\ &\bigtriangledown_bL(w,b,\xi,\alpha,\mu)=-\sum_{i=1}^N\alpha_iy_i=0\\ &\bigtriangledown_{\xi_i}L(w,b,\xi,\alpha,\mu)=C-\alpha_i-\mu_i=0 \end{aligned} $$

得到

$$ \begin{aligned} &w=\sum_{i=1}^N\alpha_iy_ix_i\\ &\sum_{i=1}^N\alpha_iy_i=0\\ &C-\alpha_i-\mu_i=0\\ \end{aligned} $$

回代到 Lagrange 函数,得到

$$ \mathop{\text{min}}_{w,b,\xi}L(w,b,\xi,\alpha,\mu)=-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(x_i\cdot x_j)+\sum_{i=1}^N\alpha_i $$

对偶问题:

$$ \begin{aligned} &\mathop{\text{max}}_{\alpha}\quad -\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(x_i\cdot x_j)+\sum_{i=1}^N\alpha_i\\ &\text{s.t.}\ \quad \sum_{i=1}^N\alpha_iy_i=0\\ &\ \\ \quad \quad C-\alpha_i-\mu_i=0\\ &\ \\ \quad \quad \alpha_i\geqslant0\\ &\ \\ \quad \quad \mu_i\geqslant0,\ i=1,2,...,N\\ \end{aligned} $$

得到如下等价问题:

$$ \begin{array}{ll} \min _{\alpha} & \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_{i} \alpha_{j} y_{i} y_{j}\left(x_{i} \cdot x_{j}\right)-\sum_{i=1}^{N} \alpha_{i} \\ \text { s.t. } & \sum_{i=1}^{N} \alpha_{i} y_{i}=0 \\ & 0 \leq \alpha_{i} \leq C, \quad i=1,2, \cdots, N \end{array} $$

对于样本向量 $\mathbf {x}_i$ ,令 $u_i = \mathbf {w} \cdot \mathbf {x}_i + b$

最优解需要满足 KKT 条件:

$$ \begin{gathered} \alpha_{i}=0 \Leftrightarrow y_{i} u_{i} \geq 1 \\ 0<\alpha_{i}<C \Leftrightarrow y_{i} u_{i}=1 \\ \alpha_{i}=C \Leftrightarrow y_{i} u_{i} \leq 1 \end{gathered} $$

若 $\exists \alpha _i$ 不满足 KKT 条件,则需要更新 $\alpha _i$ 使其满足。

为了满足条件 $ \sum_{i=1}^N\alpha_iy_i=0$ ,必须成对更新 $\alpha $,以保持等于 0.

更新的方法如下:

$$ \alpha_{1}^{\text {new }} y_{1}+\alpha_{2}^{\text {new }} y_{2}=\alpha_{1}^{\text {old }} y_{1}+\alpha_{2}^{\text {old }} y_{2}=\zeta $$

其中 $\zeta$ 为常数。$\alpha _2 ^{\text {new}}$ 的范围为 $L \leq \alpha _{2}^{\text {new}} \leq H$ 其中 $L=\max (0,-\zeta), \quad H=\min (C, C-\zeta)$

序列最小优化(SMO)算法

1996 年,John Platt 发明了 SMO 算法,也称 Platt 算法。它通过将大优化问题分解为多个小优化问题来求解。

SMO 通过求解出一系列 $\alpha $ 和 b,然后根据这些 $\alpha $ 计算权重向量 $\mathbf {w}$

辅助函数:

下面的代码摘自仓库 Jack-Cherish/Machine-Learning,但真正作者应该是 Peter Harrington

"""
读取数据
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

SMO 简化版代码如下:

"""
函数说明:简化版 SMO 算法
Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
Annotated by:
    Jack Cui
"""
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

完整 SMO 算法:

# -*-coding:utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random

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 - 包含核函数信息的元组
	Returns:
	    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 loadDataSet(fileName):
	"""
	读取数据
	Parameters:
	    fileName - 文件名
	Returns:
	    dataMat - 数据矩阵
	    labelMat - 数据标签
	"""
	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):
	"""
	计算误差
	Parameters:
		oS - 数据结构
		k - 标号为 k 的数据
	Returns:
	    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 selectJrand(i, m):
	"""
	函数说明:随机选择 alpha_j 的索引值
	Parameters:
	    i - alpha_i 的索引值
	    m - alpha 参数个数
	Returns:
	    j - alpha_j 的索引值
	"""
	j = i                                 #选择一个不等于 i 的 j
	while (j == i):
		j = int(random.uniform(0, m))
	return j

def selectJ(i, oS, Ei):
	"""
	内循环启发方式 2
	Parameters:
		i - 标号为 i 的数据的索引值
		oS - 数据结构
		Ei - 标号为 i 的数据误差
	Returns:
	    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, 并更新误差缓存
	Parameters:
		oS - 数据结构
		k - 标号为 k 的数据的索引值
	Returns:
	"""
	Ek = calcEk (oS, k)										# 计算 Ek
	oS.eCache [k] = [1,Ek]									# 更新误差缓存


def clipAlpha(aj,H,L):
	"""
	修剪 alpha_j
	Parameters:
	    aj - alpha_j 的值
	    H - alpha 上限
	    L - alpha 下限
	Returns:
	    aj - 修剪后的 alpah_j 的值
	"""
	if aj > H:
		aj = H
	if L > aj:
		aj = L
	return aj

def innerL(i, oS):
	"""
	优化的 SMO 算法
	Parameters:
		i - 标号为 i 的数据的索引值
		oS - 数据结构
	Returns:
		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)):
	"""
	完整的线性 SMO 算法
	Parameters:
		dataMatIn - 数据矩阵
		classLabels - 数据标签
		C - 松弛变量
		toler - 容错率
		maxIter - 最大迭代次数
		kTup - 包含核函数信息的元组
	Returns:
		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 testRbf(k1 = 1.3):
	"""
	测试函数
	Parameters:
		k1 - 使用高斯核函数的时候表示到达率
	Returns:
	"""
	dataArr,labelArr = loadDataSet ('testSetRBF.txt')						# 加载训练集
	b,alphas = smoP (dataArr, labelArr, 200, 0.0001, 100, ('rbf', k1))		# 根据训练集计算 b 和 alphas
	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,:],('rbf', k1))				# 计算各个点的核
		predict = kernelEval.T * np.multiply (labelSV,alphas [svInd]) + b 	# 根据支持向量的点,计算超平面,返回预测结果
		if np.sign (predict) != np.sign (labelArr [i]): errorCount += 1		# 返回数组中各元素的正负符号,用 1 和 - 1 表示,并统计错误个数
	print ("训练集错误率: %.2f%%" % ((float (errorCount)/m)*100)) 			# 打印错误率
	dataArr,labelArr = loadDataSet ('testSetRBF2.txt') 						# 加载测试集
	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,:],('rbf', k1)) 				# 计算各个点的核
		predict=kernelEval.T * np.multiply (labelSV,alphas [svInd]) + b 		# 根据支持向量的点,计算超平面,返回预测结果
		if np.sign (predict) != np.sign (labelArr [i]): errorCount += 1    	# 返回数组中各元素的正负符号,用 1 和 - 1 表示,并统计错误个数
	print ("测试集错误率: %.2f%%" % ((float (errorCount)/m)*100)) 			# 打印错误率


def showDataSet(dataMat, labelMat):
	"""
	数据可视化
	Parameters:
	    dataMat - 数据矩阵
	    labelMat - 数据标签
	Returns:
	"""
	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()

if __name__ == '__main__':
	testRbf()

本文持续更新中……

参考

(1)Support vector machines.

(2)Interactive demo of Support Vector Machines (SVM) (jgreitemann.github.io) 在线玩 SVM!