支持向量机(SVM)推导以及代码实现

SVM 是一个非常优雅的算法，具有完善的数学理论，虽然如今工业界用到的不多，但还是决定花点时间去写篇文章整理一下。

1. 支持向量

1.1 线性可分

首先我们先来了解下什么是线性可分。

在二维空间上，两类点被一条直线完全分开叫做线性可分。

严格的数学定义是：

$D_0$ 和 $D_1$ 是 n 维欧氏空间中的两个点集。如果存在 n 维向量 w 和实数 b，使得所有属于 $D_0$ 的点 $x_i$ 都有 $wx_i + b > 0$ ，而对于所有属于 $D_1$ 的点 $x_j$ 则有 $wx_j + b < 0$ ，则我们称 $D_0$ 和 $D_1$ 线性可分。

1.2 最大间隔超平面

从二维扩展到多维空间中时，将 $D_0$ 和 $D_1$ 完全正确地划分开的 $wx+b=0$ 就成了一个超平面。

为了使这个超平面更具鲁棒性，我们会去找最佳超平面，以最大间隔把两类样本分开的超平面，也称之为最大间隔超平面。

两类样本分别分割在该超平面的两侧；
两侧距离超平面最近的样本点到超平面的距离被最大化了。

1.3 支持向量

样本中距离超平面最近的一些点，这些点叫做支持向量。

1.4 SVM 最优化问题

SVM 想要的就是找到各类样本点到超平面的距离最远，也就是找到最大间隔超平面。任意超平面可以用下面这个线性方程来描述：

$w^Tx+b=0 \\$

二维空间点 $(x,y)$ 到直线 $Ax+By+C=0$ 的距离公式是：

$\frac{|Ax+By+C|}{\sqrt{A^2+B^2}} \\$

扩展到 n 维空间后，点 $x=(x_1,x_2…x_n)$ 到直线 $w^Tx+b=0$ 的距离为：

$\frac{|w^Tx+b|}{||w||} \\$

其中 $||w||=\sqrt{w_1^2+…w_n^2}$ 。

如图所示，根据支持向量的定义我们知道，支持向量到超平面的距离为 d，其他点到超平面的距离大于 d。

于是我们有这样的一个公式：

$\left\{ \begin{aligned} \frac{w^Tx+b}{||w||} &\geq d \quad y=1 \\ \frac{w^Tx+b}{||w||} &\leq -d \quad y=-1 \end{aligned} \right. \\$

稍作转化可以得到：

$\left\{ \begin{aligned} \frac{w^Tx+b}{||w||d} &\geq 1 \quad y=1 \\ \frac{w^Tx+b}{||w||d} &\leq -1 \quad y=-1 \end{aligned} \right. \\$

$||w|| d$ 是正数，我们暂且令它为 1（之所以令它等于 1，是为了方便推导和优化，且这样做对目标函数的优化没有影响），故：

$\left\{ \begin{aligned} w^Tx+b &\geq 1 \quad y=1 \\ w^Tx+b &\leq -1 \quad y=-1 \end{aligned} \right. \\$

将两个方程合并，我们可以简写为：

$y(w^Tx+b) \geq 1 \\$

至此我们就可以得到最大间隔超平面的上下两个超平面：

每个支持向量到超平面的距离可以写为：

$d=\frac{|w^Tx+b|}{||w||} \\$

由上述 $y(w^Tx+b) > 1 > 0$ 可以得到 $y(w^Tx+b) = |w^Tx+b|$ ，所以我们得到：

$d = \frac{y(w^Tx+b)}{||w||} \\$

最大化这个距离：

$\max 2* \frac{y(w^Tx+b)}{||w||} \\$

这里乘上 2 倍也是为了后面推导，对目标函数没有影响。刚刚我们得到支持向量 $y(w^Tx+b) = 1$ ，所以我们得到：

$\max \frac{2}{||w||} \\$

再做一个转换：

$\min \frac{1}{2}||w|| \\$

为了方便计算（去除 $||w||$ 的根号），我们有：

$\min \frac{1}{2}||w||^2\\$

所以得到的最优化问题是：

$\min \frac{1}{2} ||w||^2 \ s.t. \quad y_i（w^Tx_i+b）\geq 1 \\$

2. 对偶问题

2.1 拉格朗日乘数法

2.1.1 等式约束优化问题

本科高等数学学的拉格朗日程数法是等式约束优化问题：

$\min f(x_{1} ,x_{2} ,...,x_{n} ) \\ s.t. \quad h_{k} (x_{1} ,x_{2} ,...,x_{n} )=0 \quad k =1,2,...,l\\$

我们令 $L(x,\lambda ) = f(x) + \sum\limits_{k = 1}^l \lambda _k h_k (x)$ ，函数 $L(x,y)$ 称为 Lagrange 函数，参数 $\lambda$ 称为 Lagrange 乘子没有非负要求。

利用必要条件找到可能的极值点：

$\left\{ \begin{aligned} \frac{\partial L}{\partial x_i} = 0 \quad i=1,2,...,n \\ \frac{\partial L}{\partial \lambda_k} = 0 \quad k=1,2,...,l \end{aligned} \right. \\$

具体是否为极值点需根据问题本身的具体情况检验。这个方程组称为等式约束的极值必要条件。

等式约束下的 Lagrange 乘数法引入了 $l$ 个 Lagrange 乘子，我们将 $x_{i}$ 与 $\lambda_{k}$ 一视同仁，把 $\lambda_{k}$ 也看作优化变量，共有 $(n+l)$ 个优化变量。

2.1.2 不等式约束优化问题

而我们现在面对的是不等式优化问题，针对这种情况其主要思想是将不等式约束条件转变为等式约束条件，引入松弛变量，将松弛变量也是为优化变量。

以我们的例子为例：

$min f(w) = min\frac{1}{2} ||w||^2 \\ s.t. \quad g_i(w) = 1 - y_i(w^Tx_i+b)\leq 0 \\$

我们引入松弛变量 $a_i^2$ 得到 $h_i(w,a_i) = g_i(w) + a_i^2 = 0$ 。这里加平方主要为了不再引入新的约束条件，如果只引入 $a_i$ 那我们必须要保证 $a_i \geq 0$ 才能保证 $h_i(w,a_i) = 0$ ，这不符合我们的意愿。

由此我们将不等式约束转化为了等式约束，并得到 Lagrange 函数：

$\begin{aligned} L(w,\lambda,a) &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i h_i (w) \\ &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i [g_i(w) + a_i^2] \quad \lambda_i \geq 0 \end{aligned} \\$

由等式约束优化问题极值的必要条件对其求解，联立方程：

$\left\{ \begin{aligned} \frac{\partial L}{\partial w_i} &= \frac{\partial f}{\partial w_i} + \sum\limits_{i=1}^{n} \lambda_i \frac{\partial g_i}{\partial w_i}= 0, \\ \frac{\partial L}{\partial a_i} &= 2 \lambda_i a_i = 0, \\ \frac{\partial L}{\partial \lambda_i} &= g_i(w) + a_i^2 = 0, \\ \lambda_i &\geq 0 \end{aligned} \right. \\$

（为什么取 $\lambda_i \geq 0$ ，可以通过几何性质来解释，有兴趣的同学可以查下 KKT 的证明）。

针对 $\lambda_i a_i = 0$ 我们有两种情况：

情形一： $\lambda_i = 0, a_i \neq 0$

由于 $\lambda_i=0$ ，因此约束条件 $g_i(w)$ 不起作用，且 $g_i(w)<0$

情形二： $\lambda_i \neq 0, a_i = 0$

此时 $g_i(w)=0$ 且 $\lambda_i>0$ ，可以理解为约束条件 $g_i(w)$ 起作用了，且 $g_i(w)=0$

综合可得： $\lambda_ig_i(w)=0$ ，且在约束条件起作用时 $\lambda_i>0,g_i(w)=0$ ；约束不起作用时 $\lambda_i = 0,g_i(w) < 0$

由此方程组转换为：

$\left\{ \begin{aligned} \frac{\partial L}{\partial w_i} &= \frac{\partial f}{\partial w_i} + \sum\limits_{j=1}^{n} \lambda_j \frac{\partial g_j}{\partial w_i}= 0, \\ \lambda_ig_i(w) &= 0, \\ g_i(w)&\leq 0 \\ \lambda_i &\geq 0 \end{aligned} \right. \\$

以上便是不等式约束优化优化问题的 KKT(Karush-Kuhn-Tucker) 条件， $\lambda_i$ 称为 KKT 乘子。

这个式子告诉了我们什么事情呢？

直观来讲就是，支持向量 $g_i(w)=0$ ，所以 $\lambda_i > 0$ 即可。而其他向量 $g_i(w)<0, \lambda_i=0$ 。

我们原本问题时要求： $min \frac{1}{2} ||w||^2$ ，即求 $minL(w,\lambda,a)$

$\begin{aligned} L(w,\lambda,a) &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i [g_i(w) + a_i^2] \quad \\ &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i g_i(w) + \sum\limits_{i = 1}^n \lambda_i a_i^2 \end{aligned} \\$

由于 $\sum\limits_{i = 1}^n \lambda_i a_i^2 \geq 0$ ，故我们将问题转换为： $minL(w,\lambda)$ ：

$L(w,\lambda)={f(w)} + \sum\limits_{i = 1}^n \lambda_i g_i(w) \\$

假设找到了最佳参数是的目标函数取得了最小值 p。即 $\frac{1}{2} ||w||^2 =p$ 。而根据 $\lambda_{i} \geq 0$ ，可知 $\sum\limits_{i = 1}^n \lambda_i g_i(w) \leq 0$ ，因此 $L(w,\lambda) \leq p$ ，为了找到最优的参数 ${\lambda}$ ，使得 $L(w,\lambda)$ 接近 p，故问题转换为出 $\max\limits_{\lambda}L(w,\lambda)$ 。

故我们的最优化问题转换为：

$\min\limits_w \max\limits_{\lambda} L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

出了上面的理解方式，我们还可以有另一种理解方式：由于 $\lambda_i \geq 0$ ，

$\max\limits_{\lambda} L(w, \lambda) = \left\{ \begin{aligned} \infty \quad g_i(w) \geq 0 \\ \frac{1}{2} {||w||^2} \quad g_i(w) \leq 0 \end{aligned} \right. \\$

所以 $\min(\infty, \frac{1}{2} {||w||^2}) = \frac{1}{2} {||w||^2}$ ，所以转化后的式子和原来的式子也是一样的。

2.2 强对偶性

对偶问题其实就是将：

$\min\limits_w \max\limits_{\lambda} L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

变成了：

$\max\limits_{\lambda} \min\limits_w L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

假设有个函数 $f$ 我们有：

$\min\max f \geq \max\min f \\$

也就是说，最大的里面挑出来的最小的也要比最小的里面挑出来的最大的要大。这关系实际上就是弱对偶关系，而强对偶关系是当等号成立时，即：

$\min\max f = \max\min f \\$

如果 $f$ 是凸优化问题，强对偶性成立。而我们之前求的 KKT 条件是强对偶性的充要条件。

3. SVM 优化

我们已知 SVM 优化的主问题是：

$\min\limits_{w} \frac{1}{2} ||w||^2 \\ s.t.\quad g_i(w,b) = 1 - y_i（w^Tx_i+b）\leq 0, \quad i=1,2,...,n \\$

那么求解线性可分的 SVM 的步骤为：

步骤 1：

构造拉格朗日函数：

$\min\limits_{w,b}\max\limits_{\lambda} L(w,b,\lambda)= \frac{1}{2}{||w||}^2 + \sum\limits_{i = 1}^n \lambda_i [1-y_i(w^Tx_i+b)] \\ s.t. \quad \lambda_i \geq 0 \\$

步骤 2：

利用强对偶性转化：

$\max\limits_{\lambda}\min\limits_{w,b} L(w,b,\lambda) \\$

现对参数 w 和 b 求偏导数：

$\begin{aligned} \frac{\partial L}{\partial w} &= w - \sum_{i=1}^{n}\lambda_ix_iy_i = 0 \\ \frac{\partial L}{\partial b} &= \sum_{i=1}^{n}\lambda_iy_i = 0 \\ \end{aligned} \\$

得到：

$\begin{aligned} \sum_{i=1}^{n}\lambda_ix_iy_i &= w\\ \sum_{i=1}^{n}\lambda_iy_i &= 0 \\ \end{aligned} \\$

我们将这个结果带回到函数中可得：

$\begin{aligned} L(w,b,\lambda) &= \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) + \sum_{i = 1}^n \lambda_i - \sum_{i = 1}^n \lambda_i y_i(\sum_{j = 1}^n \lambda_j y_j (x_i \cdot x_j) + b) \\ &= \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) + \sum_{i = 1}^n \lambda_i - \sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j)-\sum_{i = 1}^n \lambda_i y_i b \\ &= \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \end{aligned} \\$

也就是说：

$\min\limits_{w,b}L(w,b,\lambda) = \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\$

步骤 3：

由步骤 2 得：

$\max\limits_{\lambda} [\sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j)] \\ s.t. \quad \sum_{i=1}^{n}\lambda_iy_i = 0 \quad \lambda_i \geq 0 \\$

我们可以看出来这是一个二次规划问题，问题规模正比于训练样本数，我们常用 SMO(Sequential Minimal Optimization) 算法求解。

SMO(Sequential Minimal Optimization)，序列最小优化算法，其核心思想非常简单：每次只优化一个参数，其他参数先固定住，仅求当前这个优化参数的极值。我们来看一下 SMO 算法在 SVM 中的应用。

我们刚说了 SMO 算法每次只优化一个参数，但我们的优化目标有约束条件： $\sum\limits_{i=1}^{n}\lambda_iy_i = 0$ ，没法一次只变动一个参数。所以我们选择了一次选择两个参数。具体步骤为：

选择两个需要更新的参数 $\lambda_i$ 和 $\lambda_j$ ，固定其他参数。于是我们有以下约束：

这样约束就变成了：

$\lambda_i y_i+\lambda_j y_j = c \quad \lambda_i \geq 0,\lambda_j \geq 0 \\$

其中 $c=-\sum\limits_{k \ne i,j}\lambda_ky_k$ ，由此可以得出 $\lambda_j=\frac{c-\lambda_iy_i}{y_j}$ ，也就是说我们可以用 $\lambda_i$ 的表达式代替 $\lambda_{j}$ 。这样就相当于把目标问题转化成了仅有一个约束条件的最优化问题，仅有的约束是 $\lambda_i \geq 0$ 。

2. 对于仅有一个约束条件的最优化问题，我们完全可以在 $\lambda_{i}$ 上对优化目标求偏导，令导数为零，从而求出变量值 $\lambda_{i_{new}}$ ，然后根据 $\lambda_{i_{new}}$ 求出 $\lambda_{j_{new}}$ 。

3. 多次迭代直至收敛。

通过 SMO 求得最优解 $\lambda^*$ 。

步骤 4 ：

我们求偏导数时得到：

$w = \sum_{i=1}^m \lambda_i y_i x_i \\$

由上式可求得 w。

我们知道所有 $\lambda_i > 0$ 对应的点都是支持向量，我们可以随便找个支持向量，然后带入： $y_s(wx_s+b) = 1$ ，求出 b 即可，

两边同乘 $y_s$ ，得 $y_s^2(wx_s+b) = y_s$

因为 $y_s^2=1$ ，所以： $b=y_s - wx_s$

为了更具鲁棒性，我们可以求得支持向量的均值：

$b = \frac{1}{|S|}\sum_{s \in S} (y_s -wx_s) \\$

步骤 5： w 和 b 都求出来了，我们就能构造出最大分割超平面： $w^Tx+b=0$

分类决策函数： $f(x)=sign(w^Tx+b)$

其中 $sign( \cdot )$ 为阶跃函数：

$sign(x) = \left\{ \begin{aligned} -1 \quad x<0 \\ 0 \quad x=0 \\ 1 \quad x>0 \end{aligned} \right. \\$

将新样本点导入到决策函数中既可得到样本的分类。

4. 软间隔

4.1 解决问题

在实际应用中，完全线性可分的样本是很少的，如果遇到了不能够完全线性可分的样本，我们应该怎么办？比如下面这个：

于是我们就有了软间隔，相比于硬间隔的苛刻条件，我们允许个别样本点出现在间隔带里面，比如：

我们允许部分样本点不满足约束条件：

$1-y_i(w^Tx_i + b) \leq 0 \\$

为了度量这个间隔软到何种程度，我们为每个样本引入一个松弛变量 $\xi_{i}$ ，令 $\xi_{i} \geq 0$ ，且 $1 - y_i(w^Tx_i + b)-\xi_i \leq 0$ 。对应如下图所示：

4.2 优化目标及求解

增加软间隔后我们的优化目标变成了：

$\min\limits_{w} \frac{1}{2} ||w||^2 + C\sum_{i=1}^{m}\xi_i \\ s.t.\quad g_i(w,b) = 1 - y_i(w^Tx_i+b) - \xi_i\leq 0, \quad \xi_i \geq 0, \quad i=1,2,...,n \\$

其中 C 是一个大于 0 的常数，可以理解为错误样本的惩罚程度，若 C 为无穷大， $\xi_{i}$ 必然无穷小，如此一来线性 SVM 就又变成了线性可分 SVM；当 C 为有限值的时候，才会允许部分样本不遵循约束条件。

接下来我们将针对新的优化目标求解最优化问题：

步骤 1：

构造拉格朗日函数：

$\min\limits_{w,b,\xi}\max\limits_{\lambda, \mu} L(w,b,\xi,\lambda,\mu)= \frac{1}{2}{||w||}^2 + C\sum_{i=1}^{m}\xi_i+ \sum\limits_{i = 1}^n \lambda_i [1-\xi_i-y_i(w^Tx_i+b)] - \sum_{i=1}^{n}\mu_i\xi_i \\ s.t. \quad \lambda_i \geq 0 \quad \mu_i \geq 0 \\$

其中 $\lambda_{i}$ 和 $\mu_{i}$ 是拉格朗日乘子，w、b 和 $\xi_{i}$ 是主问题参数。

根据强对偶性，将对偶问题转换为：

$\max\limits_{\lambda, \mu}\min\limits_{w,b,\xi} L(w,b,\xi,\lambda,\mu) \\$

步骤 2：

分别对主问题参数w、b 和 $\xi_{i}$ 求偏导数，并令偏导数为 0，得出如下关系：

$w = \sum_{i=1}^{m}\lambda_i y_i x_i \\ 0 = \sum_{i=1}^{m}\lambda_i y_i \\ C = \lambda_i + \mu_i \\$

将这些关系带入拉格朗日函数中，得到：

$\min\limits_{w,b,\xi}L(w,b,\xi,\lambda,\mu) = \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\$

最小化结果只有 $\lambda$ 而没有 $\mu$ ，所以现在只需要最大化 $\lambda$ 就好：

我们可以看到这个和硬间隔的一样，只是多了个约束条件。

然后我们利用 SMO 算法求解得到拉格朗日乘子 $\lambda^{*}$ 。

步骤 3 ：

$w = \sum_{i=1}^m \lambda_i y_i x_i \\ b = \frac{1}{|S|}\sum_{s \in S} (y_s -wx_s) \\$

然后我们通过上面两个式子求出 w 和 b，最终求得超平面 $w^Tx+b=0$ ，

这边要注意一个问题，在间隔内的那部分样本点是不是支持向量？

我们可以由求参数 w 的那个式子可看出，只要 $\lambda_{i} > 0$ 的点都能够影响我们的超平面，因此都是支持向量。

5. 核函数

5.1 线性不可分

我们刚刚讨论的硬间隔和软间隔都是在说样本的完全线性可分或者大部分样本点的线性可分。

但我们可能会碰到的一种情况是样本点不是线性可分的，比如：

这种情况的解决方法就是：将二维线性不可分样本映射到高维空间中，让样本点在高维空间线性可分，比如：

对于在有限维度向量空间中线性不可分的样本，我们将其映射到更高维度的向量空间里，再通过间隔最大化的方式，学习得到支持向量机，就是非线性 SVM。

我们用 x 表示原来的样本点，用 $\phi(x)$ 表示 x 映射到特征新的特征空间后到新向量。那么分割超平面可以表示为： $f(x)=w \phi(x)+b$ 。

对于非线性 SVM 的对偶问题就变成了：

$\min\limits_{\lambda} [\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (\phi(x_i) \cdot \phi(x_j))-\sum_{j=1}^{n}\lambda_i] \\ s.t. \quad \sum_{i=1}^{n}\lambda_iy_i = 0, \quad \lambda_i \geq 0, \quad C-\lambda_i-\mu_i=0 \\$

可以看到与线性 SVM 唯一的不同就是：之前的 $(x_i \cdot x_j)$ 变成了 $(\phi(x_i) \cdot \phi(x_j))$ 。

5.2 核函数的作用

我们不禁有个疑问：只是做个内积运算，为什么要有核函数的呢？

这是因为低维空间映射到高维空间后维度可能会很大，如果将全部样本的点乘全部计算好，这样的计算量太大了。

但如果我们有这样的一核函数 $k(x,y) = (\phi(x),\phi(y))$ ， $x_i$ 与 $x_j$ 在特征空间的内积等于它们在原始样本空间中通过函数 $k( x, y)$ 计算的结果，我们就不需要计算高维甚至无穷维空间的内积了。

举个例子：假设我们有一个多项式核函数：

$k(x,y)=(x \cdot y + 1)^2 \\$

带进样本点的后：

$k(x,y) = (\sum_{i=1}^n(x_i \cdot y_i) + 1)^2 \\$

而它的展开项是：

$\sum_{i=1}^nx_i^2y_i^2+\sum_{i=2}^n\sum_{j=1}^{i-1}(\sqrt2x_ix_j)(\sqrt2y_iy_j)+\sum_{i=1}{n}(\sqrt2x_i)(\sqrt2y_i)+1 \\$

如果没有核函数，我们则需要把向量映射成：

$x^{'} = (x_1^2,...,x_n^2,...\sqrt2x_1,...,\sqrt2x_n,1) \\$

然后在进行内积计算，才能与多项式核函数达到相同的效果。

可见核函数的引入一方面减少了我们计算量，另一方面也减少了我们存储数据的内存使用量。

5.3 常见核函数

我们常用核函数有：

线性核函数

$k(x_i,x_j) = x_i^Tx_j \\$

多项式核函数

$k(x_i,x_j) = (x_i^Tx_j)^d\\$

高斯核函数

$k(x_i,x_j) = exp(-\frac{||x_i-x_j||}{2\delta^2}) \\$

这三个常用的核函数中只有高斯核函数是需要调参的。

6. 优缺点

6.1 优点

有严格的数学理论支持，可解释性强，不依靠统计方法，从而简化了通常的分类和回归问题；
能找出对任务至关重要的关键样本（即：支持向量）；
采用核技巧之后，可以处理非线性分类/回归任务；
最终决策函数只由少数的支持向量所确定，计算的复杂性取决于支持向量的数目，而不是样本空间的维数，这在某种意义上避免了“维数灾难”。

6.2 缺点

训练时间长。当采用 SMO 算法时，由于每次都需要挑选一对参数，因此时间复杂度为 $O(N^2)$ ，其中 N 为训练样本的数量；
当采用核技巧时，如果需要存储核矩阵，则空间复杂度为 $O(N^2)$ ；
模型预测时，预测时间与支持向量的个数成正比。当支持向量的数量较大时，预测计算复杂度较高。

因此支持向量机目前只适合小批量样本的任务，无法适应百万甚至上亿样本的任务。

import random
from numpy import *


# SMO算法相关辅助中的辅助函数
# 1 解析文本数据函数，提取每个样本的特征组成向量，添加到数据矩阵
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    with open(fileName) as f:
        for line in f.readlines():
            lineArr = line.strip().split('\t')
            dataMat.append([float(lineArr[0]), float(lineArr[1])])
            labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


# 2 在样本集中采取随机选择的方法选取第二个不等于第一个alphai的
# 优化向量alphaj
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


# 3 约束范围L<=alphaj<=H内的更新后的alphaj值
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# @dataMat    ：数据列表
# @classLabels：标签列表
# @C          ：权衡因子（增加松弛因子而在目标优化函数中引入了惩罚项）
# @toler      ：容错率
# @maxIter    ：最大迭代次数
def smoSimple(dataMat, classLabels, C, toler, maxIter):
    # 将列表形式转为矩阵或向量形式
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    # 初始化b=0，获取矩阵行列
    b = 0
    m, n = shape(dataMatrix)
    # 新建一个m行1列的向量
    alphas = mat(zeros((m, 1)))
    # 迭代次数为0
    iter = 0
    while (iter < maxIter):
        # 改变的alpha对数
        alphaPairsChanged = 0
        # 遍历样本集中样本
        for i in range(m):
            # 计算支持向量机算法的预测值
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            # 计算预测值与实际值的误差
            Ei = fXi - float(labelMat[i])
            # 如果不满足KKT条件，即labelMat[i]*fXi<1(labelMat[i]*fXi-1<-toler)
            # and alpha<C 或者labelMat[i]*fXi>1(labelMat[i]*fXi-1>toler)and alpha>0
            if ((labelMat[i] * Ei < -toler) and (alpha < C)) or ((labelMat[i] * Ei > toler) and (alpha[i] > 0)):
                # 随机选择第二个变量alphaj
                j = selectJrand(i, m)
                # 计算第二个变量对应数据的预测值
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :]).T) + b
                # 计算与测试与实际值的差值
                Ej = fXj - float(label[j])
                # 记录alphai和alphaj的原始值，便于后续的比较
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # 如何两个alpha对应样本的标签不相同
                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
                # 根据公式计算未经剪辑的alphaj
                # ------------------------------------------
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - ataMatrix[i, :] * dataMatrix[i, :].T - dataMatrix[j,:] * dataMatrix[j, :].T
                # 如果eta>=0,跳出本次循环
                if eta >= 0: print("eta>=0");continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                # ------------------------------------------
                # 如果改变后的alphaj值变化不大，跳出本次循环
                if (abs(alphas[j] - alphaJold) < 0.00001): print("j not moving enough");continue
                # 否则，计算相应的alphai值
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                # 再分别计算两个alpha情况下对于的b值
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMat[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
                # 如果0<alphai<C,那么b=b1
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                # 否则如果0<alphai<C,那么b=b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                # 否则，alphai，alphaj=0或C
                else:
                    b = (b1 + b2) / 2.0
                # 如果走到此步，表面改变了一对alpha值
                alphaPairsChanged += 1
                print("iter: %d i:%d,paird changed %d" % (iter, i, alphaPairsChanged))
        # 最后判断是否有改变的alpha对，没有就进行下一次迭代
        if (alphaPairsChanged == 0):
            iter += 1
        # 否则，迭代次数置0，继续循环
        else:
            iter = 0
        print("iteration number: %d" % iter)
    # 返回最后的b值和alpha向量
    return b, alphas

#启发式SMO算法的支持函数
#新建一个类的收据结构，保存当前重要的值
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))


# 格式化计算误差的函数，方便多次调用
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek


# 修改选择第二个变量alphaj的方法
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    # 将误差矩阵每一行第一列置1，以此确定出误差不为0
    # 的样本
    oS.eCache[i] = [1, Ei]
    # 获取缓存中Ei不为0的样本对应的alpha列表
    validEcacheList = nonzero(oS.Cache[:, 0].A)[0]
    # 在误差不为0的列表中找出使abs(Ei-Ej)最大的alphaj
    if (len(validEcacheList) > 0):
        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:
        # 否则，就从样本集中随机选取alphaj
        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]


        # 内循环寻找alphaj

    def innerL(i, oS):
        # 计算误差
        Ei = calcEk(oS, i)
        # 违背kkt条件
        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)
            alphaIold = alphas[i].copy()
            alphaJold = 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.alphas[j] + oS.alphas[i])
            if L == H: print("L==H");return 0
            # 计算两个alpha值
            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] - alphaJold) < 0.00001):
                print("j not moving enough")
                return 0
            oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
            updateEk(oS, i)
            # 在这两个alpha值情况下，计算对应的b值
            # 注，非线性可分情况，将所有内积项替换为核函数K[i,j]
            b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
            b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * 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:
                oS.b = (b1 + b2) / 2.0
            # 如果有alpha对更新
            return 1
            # 否则返回0
        else: return 0

# SMO外循环代码
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    # 保存关键数据
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # 选取第一个变量alpha的三种情况，从间隔边界上选取或者整个数据集
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 没有alpha更新对
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print("fullSet,iter: %d i:%d,pairs changed %d" %(iter, i, alphaPairsChanged))
            iter += 1
        else:
            # 统计alphas向量中满足0<alpha<C的alpha列表
            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,pairs changed %d"% (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
            # 如果本次循环没有改变的alpha对，将entireSet置为true，
            # 下个循环仍遍历数据集
        elif (alphaPairsChanged == 0):
            entireSet = True
            print("iteration number: %d" % iter)
    return oS.b, oS.alphas

# 求出了alpha值和对应的b值，就可以求出对应的w值，以及分类函数值
def predict(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels)
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    result = dataArr[0] * mat(ws) + b
    return sign(result)