机器学习之支持向量机SVM

函数间隔和几何间隔

一般来讲，一个样本点距离分隔超平面的远近可以表示分类预测的确信程度。为了能够示分类预测的确信程度，我们分别定义函数间隔（Functional Margin）和几何间隔（Geometric Margin）

1.函数间隔

对于给定的训练数据集$(X^{(1)},y^{(1)}),(X^{(2)},y^{(2)}),...,(X^{(m)},y^{(m)})$(X(1),y(1)),(X(2),y(2)),...,(X(m),y(m))和分隔超平面
$W \cdot X^{(i)}+b$W⋅X(i)+b，定义分超平面关于样本点$(X^{(i)},y^{(i)})$(X(i),y(i))的函数间隔为：
$\hat{\gamma}_{i}=y^{(i)}\left(W \cdot X^{(i)}+b\right)$γ^​i​=y(i)(W⋅X(i)+b)
说明：当预测值$W \cdot X^{(i)}+b$W⋅X(i)+b与$y^{(i)}$y(i)同号时，表明分类是正确的；$\hat{\gamma}_{i}$γ^​i​表示距离分割平面的远近，越远确信度越高。

函数间隔可以表示分类预测的正确性和确定性。但是，在分隔超平面中，如果其参数W 和b同时扩大为原来的2倍，这对分隔超平面来说，并没有任何改变，但是对于函数间隔 $\hat{\gamma}_{i}$γ^​i​来说，即扩大为原来的2倍。为了解决这样的问题，我们引入何间隔。

2.几何间隔

对于给定的训练数据集$(X^{(1)},y^{(1)}),(X^{(2)},y^{(2)}),...,(X^{(m)},y^{(m)})$(X(1),y(1)),(X(2),y(2)),...,(X(m),y(m))和分隔超平面
$W \cdot X^{(i)}+b$W⋅X(i)+b，定义分超平面关于样本点$(X^{(i)},y^{(i)})$(X(i),y(i))的几何间隔为：
$\gamma_{i}=y^{(i)}\left(\frac{W}{\|W\|} \cdot X^{(i)}+\frac{b}{\|W\|}\right)$γi​=y(i)(∥W∥W​⋅X(i)+∥W∥b​)
从上面的定义不难发现，几何间隔其实就是样本到分隔超平面的距离。对于几何间隔和函数间隔，有如下的关系:
$\gamma_{i}=\frac{\hat{\gamma}_{i}}{\|W\|}$γi​=∥W∥γ^​i​​

支持向量机SVM

在支持向量机（Support Vector Machines，SVM）中，求解出的分隔超平面不仅能够正确划分训练数集，而且要求几何间隔最大。

SVM模型的建立

1.硬间隔模型

$\max _{W, b} \gamma$W,bmax​γ

$y^{(i)}\left(\frac{W}{\|W\|} \cdot X^{(i)}+\frac{b}{\|W\|}\right) \geqslant \gamma, \quad i=1,2, \cdots, m$y(i)(∥W∥W​⋅X(i)+∥W∥b​)⩾γ,i=1,2,⋯,m

利用几何间隔与函数间隔的关系，原问题可以转化为：
$\max _{W, b} \frac{\hat{\gamma}}{\|W\|}$W,bmax​∥W∥γ^​​

$y^{(i)}\left(W \cdot X^{(i)}+b\right) \geqslant \hat{\gamma}, \quad i=1,2, \cdots, m$y(i)(W⋅X(i)+b)⩾γ^​,i=1,2,⋯,m

在函数间隔中，函数间隔$\hat{\gamma}_{i}$γ^​i​ 的取值并不影响到最优问题的解，如上所述，当参数W 和b同时扩大为原来的2倍，函数间隔 $\hat{\gamma}_{i}$γ^​i​也会同时扩大为原来的2倍，这对于上述的优化问题和约束条件并没有影响，因此，可以取$\hat{\gamma}_{i}$γ^​i​=1(便于运算)，故原问题可转化为：
$\max _{W, b} \frac{1}{\|W\|}$W,bmax​∥W∥1​

$y^{(i)}\left(W \cdot X^{(i)}+b\right) \geqslant1, \quad i=1,2, \cdots, m$y(i)(W⋅X(i)+b)⩾1,i=1,2,⋯,m

再将求解最大化问题转换为求解最小化问题，则上述的优化问题转换成：
$\min _{W, b} \frac{1}{2}\|W\|^{2}$W,bmin​21​∥W∥2

$\text{s.t. } y^{(i)}\left(W \cdot X^{(i)}+b\right)-1 \geqslant 0, \quad i=1,2, \cdots, m$s.t. y(i)(W⋅X(i)+b)−1⩾0,i=1,2,⋯,m

2.软间隔模型（松弛变量&惩罚因子）

对于一个数据集，可能其中存在部分的特异点(噪声)，但是将些特异点除去后，剩下的大部分的样本点组成的集合是线性可分的。
对于线性不可分的某些样本$(X^{(i)},y^{(i)})$(X(i),y(i))意味着其不能满足函数间隔大于或等于1 的约束条件，为了解决这个题，可以对每个样本$(X^{(i)},y^{(i)})$(X(i),y(i))引进一个松弛变量$\xi_{i}$ξi​ ≥0，使得函数间隔加上松弛变量大于或等于1，这样，约束条变为：
$y^{(i)}\left(W \cdot X^{(i)}+b\right) \geqslant 1-\xi_{i}$y(i)(W⋅X(i)+b)⩾1−ξi​

同时，对每个松弛变量$\xi_{i}$ξi​，支付一个代价C，此时，目标函数变为：
$\frac{1}{2}\|W\|^{2}+C \sum_{i=1}^{m} \xi_{i}$21​∥W∥2+Ci=1∑m​ξi​

此时优化目标变为：
$\min _{W, b, \xi} \frac{1}{2}\|W\|^{2}+C \sum_{i=1}^{m} \xi_{i}$W,b,ξmin​21​∥W∥2+Ci=1∑m​ξi​

$\text{s.t.}\quad y_{i}\left(W \cdot X^{(i)}+b\right) \geqslant 1-\xi_{i}, \quad i=1,2, \cdots, m \\ \xi_{i} \geq 0, \quad i=1,2, \cdots, m$s.t.yi​(W⋅X(i)+b)⩾1−ξi​,i=1,2,⋯,mξi​≥0,i=1,2,⋯,m

松弛变量

• 并非所有的样本点都有一个松弛变量与其对应
• 松弛变量的值实际上标示出了对应的点到底离群有多远，值越大，点就越远

惩罚因子

• 惩罚因子C决定了你对离群点的重视程度，显然当所有离群点的松弛变量的和一定时，你定的C越大，对目标函数的损失也越大，此时就暗示着你非常不愿意放弃这些离群点，最极端的情况是你把C定为无限大，这样只要稍有一个点离群，目标函数的值马上变成无限大，马上让问题变成无解，这就退化成了之前的硬间隔问题
• 惩罚因子C不是一个变量，整个优化问题在解的时候，C是一个你必须事先指定的值，指定这个值以后，解一下，得到一个分类器，然后用测试数据看看结果怎么样，如果不够好，换一个C的值，再解一次优化问题，得到另一个分类器，再看看效果，如此就是一个参数调优的过程，但这和优化问题本身决不是一回事，优化问题在解的过程中，C一直是定值
• 数据集偏斜（unbalanced）【它指的是参与分类的两个类别（也可以指多个类别）样本数量差异很大】要求我们需要给样本数量少的类别更大的惩罚因子，表示我们重视这类样本

模型的求解

1.拉格朗日对偶性

$\min _{W, b, \xi} \frac{1}{2}\|W\|^{2}+C \sum_{i=1}^{m} \xi_{i}$W,b,ξmin​21​∥W∥2+Ci=1∑m​ξi​

$\text{s.t.}\quad y_{i}\left(W \cdot X^{(i)}+b\right) \geqslant 1-\xi_{i}, \quad i=1,2, \cdots, m \\ \xi_{i} \geq 0, \quad i=1,2, \cdots, m$s.t.yi​(W⋅X(i)+b)⩾1−ξi​,i=1,2,⋯,mξi​≥0,i=1,2,⋯,m
对于带约束的优化问题的求解，可以使用拉格朗日乘数法，将其转化为无约束优化问题的求解。对于上述的带约束的优问题，可以转换成如下的拉格朗日函数：
${L}(w, b, \xi, \alpha, \beta)=\frac{1}{2}\|W\|^{2}+C \sum_{i=1}^{n} \xi_{i}-\sum_{i=1}^{n} \alpha_{i}\left(y^{(i)}\left(W^{T} X^{(i)}+b\right)-1+\xi_{i}\right)-\sum_{i=1}^{n} \beta_{i} \xi_{i}$L(w,b,ξ,α,β)=21​∥W∥2+Ci=1∑n​ξi​−i=1∑n​αi​(y(i)(WTX(i)+b)−1+ξi​)−i=1∑n​βi​ξi​

上述的最小优化问题即为：
$\min _{W, b, \xi} \max _{\alpha, \beta} L(W, b, \xi, \alpha, \beta)$W,b,ξmin​α,βmax​L(W,b,ξ,α,β)

根据拉格朗日对偶性，原始问题的对偶问题为：
$\max _{\alpha, \beta} \min _{W, b, \xi} L(W, b, \xi, \alpha, \beta)$α,βmax​W,b,ξmin​L(W,b,ξ,α,β)

先求$\min _{W, b, \xi} L(W, b, \xi, \alpha, \beta)$W,b,ξmin​L(W,b,ξ,α,β)

分别对$W,b,\xi$W,b,ξ求偏导
$\frac{\partial L(W, b, \xi, \alpha, \beta)}{\partial W}=W-\sum_{i=1}^{m} \alpha_{i} y^{(i)} X^{(i)}=0$∂W∂L(W,b,ξ,α,β)​=W−i=1∑m​αi​y(i)X(i)=0

$\frac{\partial L(W, b, \xi, \alpha, \beta)}{\partial b}=-\sum_{i=1}^{m} \alpha_{i} y^{(i)}=0$∂b∂L(W,b,ξ,α,β)​=−i=1∑m​αi​y(i)=0

$\frac{\partial L(W, b, \xi, \alpha, \beta)}{\partial \xi_{i}}=C-\alpha_{i}-\beta_{i}=0$∂ξi​∂L(W,b,ξ,α,β)​=C−αi​−βi​=0

代入原式得：
$\min _{W, b, \xi} L(W, b, \xi, \alpha, \beta)=-\frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \alpha_{i} \alpha_{j} y^{(i)} y^{(j)}\left(X^{(i)} \cdot X^{(j)}\right)+\sum_{i=1}^{m} \alpha_{i}$W,b,ξmin​L(W,b,ξ,α,β)=−21​i=1∑m​j=1∑m​αi​αj​y(i)y(j)(X(i)⋅X(j))+i=1∑m​αi​

即得对偶问题：
$\max _{\alpha}-\frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \alpha_{i} \alpha_{j} y^{(i)} y^{(j)}\left(X^{(i)} \cdot X^{(j)}\right)+\sum_{i=1}^{m} \alpha_{i}$αmax​−21​i=1∑m​j=1∑m​αi​αj​y(i)y(j)(X(i)⋅X(j))+i=1∑m​αi​

$\text{s.t.}\sum_{i=1}^{m} \alpha_{i} y^{(i)}=0\\C-\alpha_{i}-\beta_{i}=0\\\alpha_i \geq0\\\beta_i \geq0,i=1,2\cdots,m$s.t.i=1∑m​αi​y(i)=0C−αi​−βi​=0αi​≥0βi​≥0,i=1,2⋯,m

由于$\beta_i \geq0$βi​≥0可知$0\leq\alpha_i\leq C$0≤αi​≤C。同时将最大化问题转化为最小化问题
$\min _{\alpha} \frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \alpha_{i} \alpha_{j} y^{(i)} y^{(j)}\left(X^{(i)} \cdot X^{(j)}\right)-\sum_{i=1}^{m} \alpha_{i}$αmin​21​i=1∑m​j=1∑m​αi​αj​y(i)y(j)(X(i)⋅X(j))−i=1∑m​αi​

$\text{s.t.}\sum_{i=1}^{m} \alpha_{i} y^{(i)}=0\\0\leq\alpha_i\leq C,i=1,2,\cdots,m$s.t.i=1∑m​αi​y(i)=00≤αi​≤C,i=1,2,⋯,m

对于一个非线性可分得问题，我们可以采用核函数得方式将其转化为线性问题
在这里，我们使用高斯核函数
$K\left(X^{(i)}, X^{(j)}\right)=\exp \left(-\frac{\left\|X^{(i)}-X^{(j)}\right\|}{2 \sigma^{2}}\right)$K(X(i),X(j))=exp(−2σ2∥∥​X(i)−X(j)∥∥​​)

这样，对于非线性支持向量机得优化目标为：
$\min _{\alpha} \frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \alpha_{i} \alpha_{j} y^{(i)} y^{(j)}K\left(X^{(i)} \cdot X^{(j)}\right)-\sum_{i=1}^{m} \alpha_{i}$αmin​21​i=1∑m​j=1∑m​αi​αj​y(i)y(j)K(X(i)⋅X(j))−i=1∑m​αi​

$\text{s.t.}\sum_{i=1}^{m} \alpha_{i} y^{(i)}=0\\0\leq\alpha_i\leq C,i=1,2,\cdots,m$s.t.i=1∑m​αi​y(i)=00≤αi​≤C,i=1,2,⋯,m

终于经过一步步努力，我们已经把一个分类问题转化成一个二次规划问题了：)

2.序列最小最优化算法SMO

SMO的思想是将一个大的问题划分成一系列小的问题，通过对这些子问题的求解，达到对对偶问题的求解过程。在SMO算法中，不断将对偶问题的二次规划问题分解为只有两个变量的二次规划子问题，并对子问题进行求解。

2.1 问题转化

在此，我们先求解两个变量$\alpha_1$α1​ 和$\alpha_2$α2​的更新方法。将设取得的两个变量为$\alpha_1$α1​和$\alpha_2$α2​，其他变量为固定值，此时的优化问题为:
$\min _{\alpha_{1}, \alpha_{2}} \frac{1}{2} K_{11} \alpha_{1}^{2}+\frac{1}{2} K_{22} \alpha_{2}^{2}+y^{(1)} y^{(2)} K_{12} \alpha_{1} \alpha_{2}-\left(\alpha_{1}+\alpha_{2}\right)+y^{(1)} \alpha_{1} \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 1}+y^{(2)} \alpha_{2} \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 2}+M_{1}$α1​,α2​min​21​K11​α12​+21​K22​α22​+y(1)y(2)K12​α1​α2​−(α1​+α2​)+y(1)α1​k=3∑m​y(k)αk​Kk1​+y(2)α2​k=3∑m​y(k)αk​Kk2​+M1​

$\text{s.t.}\alpha_{1} y^{(1)}+\alpha_{2} y^{(2)}=-\sum_{k=3}^{m} \alpha_{k} y^{(k)}=M_{2}\\0\leq\alpha_k\leq C,k=1,2$s.t.α1​y(1)+α2​y(2)=−k=3∑m​αk​y(k)=M2​0≤αk​≤C,k=1,2

其中$K_{i,j}=K\left(X^{(i)}, X^{(j)}\right),M_1和M_2$Ki,j​=K(X(i),X(j)),M1​和M2​表示与$\alpha_1和\alpha+2$α1​和α+2
无关的常数项

2.2 转化为一个二元函数

$W(\alpha_1,\alpha_2)= \frac{1}{2} K_{11} \alpha_{1}^{2}+\frac{1}{2} K_{22} \alpha_{2}^{2}+y^{(1)} y^{(2)} K_{12} \alpha_{1} \alpha_{2}-\left(\alpha_{1}+\alpha_{2}\right)+y^{(1)} \alpha_{1} \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 1}+y^{(2)} \alpha_{2} \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 2}+M_{1}$W(α1​,α2​)=21​K11​α12​+21​K22​α22​+y(1)y(2)K12​α1​α2​−(α1​+α2​)+y(1)α1​k=3∑m​y(k)αk​Kk1​+y(2)α2​k=3∑m​y(k)αk​Kk2​+M1​

2.3 转化为一个一元函数

由$\alpha_{1} y^{(1)}+\alpha_{2} y^{(2)}=-\sum_{k=3}^{m} \alpha_{k} y^{(k)}=M_{2}$α1​y(1)+α2​y(2)=−∑k=3m​αk​y(k)=M2​，可以消掉一个元，从而转化为：
$W(\alpha_2)= \frac{1}{2} K_{11} {\left(M_{2}-\alpha_{2} y^{(2)}\right) }^{2}+\frac{1}{2} K_{22} \alpha_{2}^{2}+y^{(1)} y^{(2)} K_{12}{\left(M_{2}-\alpha_{2} y^{(2)}\right) } \alpha_{2}-\left({\left(M_{2}-\alpha_{2} y^{(2)}\right) }+\alpha_{2}\right)+y^{(1)} {\left(M_{2}-\alpha_{2} y^{(2)}\right) } \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 1}+y^{(2)} \alpha_{2} \sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 2}+M_{1}$W(α2​)=21​K11​(M2​−α2​y(2))2+21​K22​α22​+y(1)y(2)K12​(M2​−α2​y(2))α2​−((M2​−α2​y(2))+α2​)+y(1)(M2​−α2​y(2))k=3∑m​y(k)αk​Kk1​+y(2)α2​k=3∑m​y(k)αk​Kk2​+M1​

2.4 对一元函数求极值点

$\frac{\partial}{\partial \alpha_{2}} W\left(\alpha_{2}\right)=0$∂α2​∂​W(α2​)=0

得到：
$\begin{array}{l}{\left(K_{11}+K_{22} \alpha_{2}-2 K_{12}\right) \alpha_{2}} \\ {=\left(K_{11} M_{2}-K_{12} M_{2}-y^{(1)}+y^{(2)}+\sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 1}-\sum_{k=3}^{m} y^{(k)} \alpha_{k} K_{k 2}\right) y^{(2)}}\end{array}$(K11​+K22​α2​−2K12​)α2​=(K11​M2​−K12​M2​−y(1)+y(2)+∑k=3m​y(k)αk​Kk1​−∑k=3m​y(k)αk​Kk2​)y(2)​

令$g(X)=\sum_{i=1}^{m} y^{(i)} \alpha_{i} K\left(X^{(i)}, X\right)+b$g(X)=i=1∑m​y(i)αi​K(X(i),X)+b

$E_i=g(X^{(i)})-y^{(i)}$Ei​=g(X(i))−y(i)
假设第$i$i代，$\alpha_1^{(i)}y^{(1)}+\alpha_2^{(i)}y^{(2)}=M_2$α1(i)​y(1)+α2(i)​y(2)=M2​，则对于第$i+1$i+1代时：
$\alpha_{2}^{(i+1)}=\alpha_{2}^{(i)}+\frac{y^{(2)}\left(E_{1}-E_{2}\right)}{\eta}$α2(i+1)​=α2(i)​+ηy(2)(E1​−E2​)​

其中$\eta=K_{11}+K_{22}-K_{12}$η=K11​+K22​−K12​

2.5 解的约束条件

• $0 \leq \alpha_{i=1,2} \leq C$0≤αi=1,2​≤C
• $\alpha_1y^{(1)}+\alpha_2y^{(2)}=M_2$α1​y(1)+α2​y(2)=M2​

最优解必须要在方框内且在直线上取得,所以$L\leq\alpha_2^{(i+1)}\leq H$L≤α2(i+1)​≤H

• $当y^{(1)} 不等于 y^{(2)}$当y(1)不等于y(2)
  $L=max(0,\alpha_2^{(i)}-\alpha_1^{(i)}),H=min(C,C+\alpha_2^{(i)}-\alpha_1^{(i)})$L=max(0,α2(i)​−α1(i)​),H=min(C,C+α2(i)​−α1(i)​)
• $当y^{(1)} 等于 y^{(2)}$当y(1)等于y(2)
  $L=max(0,\alpha_2^{(i)}+\alpha_1^{(i)}-C),H=min(C,\alpha_2^{(i)}+\alpha_1^{(i)})$L=max(0,α2(i)​+α1(i)​−C),H=min(C,α2(i)​+α1(i)​)

所以，$\alpha_2^{(i+1)}$α2(i+1)​的值为：
$\alpha_{2}^{(i+1)}=\left\{\begin{array}{cc}{H} &amp; {\alpha_{2}^{(i+1)}&gt;H} \\ {\alpha_{2}^{(i+1)}} &amp; {L \leq \alpha_{2}^{(i+1)} \leq H} \\ {L} &amp; {\alpha_{2}^{(i+1)}&lt;L}\end{array}\right.$α2(i+1)​=⎩⎪⎨⎪⎧​Hα2(i+1)​L​α2(i+1)​>HL≤α2(i+1)​≤Hα2(i+1)​<L​

2.6 求解另一个变量

由$\alpha_{1}^{(i+1)} y^{(1)}+\alpha_{2}^{(i+1)} y^{(2)}=\alpha_{1}^{(i)} y^{(1)}+\alpha_{2}^{(i)} y^{(2)}=M_{2}$α1(i+1)​y(1)+α2(i+1)​y(2)=α1(i)​y(1)+α2(i)​y(2)=M2​,可推出
$\alpha_{1}^{(i+1)}=\alpha_{1}^{(i)}+y^{(1)} y^{(2)}\left(\alpha_{2}^{(i)}-\alpha_{2}^{(i+1)}\right)$α1(i+1)​=α1(i)​+y(1)y(2)(α2(i)​−α2(i+1)​)

2.7 如何选择两个变

当所有变量都满足KKT条件时，那么便是最优化问题的解.所以对于第一个变量的选择，可以选择不满足KKT条件的变量
KKT条件：
$\begin{array}{l}{\alpha_{1}=0 \Leftrightarrow y^{(1)} g\left(X^{(1)}\right) \geq 1} \\ {0&lt;\alpha_{1}&lt;C \Leftrightarrow y^{(1)} g\left(X^{(1)}\right)=1} \\ {\alpha_{1}=C \Leftrightarrow y^{(1)} g\left(X^{(1)}\right) \leq 1}\end{array}$α1​=0⇔y(1)g(X(1))≥10<α1​<C⇔y(1)g(X(1))=1α1​=C⇔y(1)g(X(1))≤1​
不满足KKT条件如下：

• 如果$y^{(1)}E_1&lt;0$y(1)E1​<0,此时，若$\alpha_1&lt;C$α1​<C，则违反KKT条件
• 如果$y^{(1)}E_1&gt;0$y(1)E1​>0,此时，若$\alpha_1&gt;0$α1​>0，则违反KKT条件

通过检查每一个样本点是否满足KKT条件，选择出第一个变量。

---

第二个变量的选择原则是要使得$\alpha_2$α2​能够发生足够大的变化，由其更新公式$\alpha_{2}^{(i+1)}=\alpha_{2}^{(i)}+\frac{y^{(2)}\left(E_{1}-E_{2}\right)}{\eta}$α2(i+1)​=α2(i)​+ηy(2)(E1​−E2​)​

由于$\alpha_1,\alpha_2$α1​,α2​是依赖于$E_1−E_2$E1​−E2​,当$E_1$E1​为正时，那么选择最小的$E_i$Ei​作为$E_2$E2​,反之，选择最大$E_i$Ei​作为$E_2$E2​

2.8 阈值b的更新与计算

当更新完成后$\alpha_1,\alpha_2$α1​,α2​，需要重新计算b
$\sum_{j=1}^{m} \alpha_{j} y^{(j)} K\left(X^{(j)}, X^{(1)}\right)+b=y^{(1)}$j=1∑m​αj​y(j)K(X(j),X(1))+b=y(1)

$b_1^{(i+1)}=y^{(1)}-\sum_{j=3}^{m}\alpha_jy^{(j)}K(X^{(j)},X^{(1)})-\alpha_1^{(i+1)}y^{(1)}K(X^{(1)},X^{(1)})-\alpha_2^{(i+1)}y^{(2)}K(X^{(2)},X^{(1)})$b1(i+1)​=y(1)−j=3∑m​αj​y(j)K(X(j),X(1))−α1(i+1)​y(1)K(X(1),X(1))−α2(i+1)​y(2)K(X(2),X(1))

SVM算法实践

import numpy as np
import pickle

class SVM:
    def __init__(self, dataSet, labels, C, toler, kernel_option):
        self.train_x = dataSet # 训练特征
        self.train_y = labels  # 训练标签
        self.C = C # 惩罚参数
        self.toler = toler     # 迭代的终止条件之一
        self.n_samples = np.shape(dataSet)[0] # 训练样本的个数
        self.alphas = np.mat(np.zeros((self.n_samples, 1))) # 拉格朗日乘子
        self.b = 0
        self.error_tmp = np.mat(np.zeros((self.n_samples, 2))) # 保存E的缓存
        self.kernel_opt = kernel_option # 选用的核函数及其参数
        self.kernel_mat = calc_kernel(self.train_x, self.kernel_opt) # 核函数的输出

def cal_kernel_value(train_x, train_x_i, kernel_option):
    '''样本之间的核函数的值
    input:  train_x(mat):训练样本
            train_x_i(mat):第i个训练样本
            kernel_option(tuple):核函数的类型以及参数
    output: kernel_value(mat):样本之间的核函数的值
            
    '''
    kernel_type = kernel_option[0] # 核函数的类型，分为rbf和其他
    m = np.shape(train_x)[0] # 样本的个数
    
    kernel_value = np.mat(np.zeros((m, 1)))
    
    if kernel_type == 'rbf': # rbf核函数
        sigma = kernel_option[1]
        if sigma == 0:
            sigma = 1.0
        for i in range(m):
            diff = train_x[i, :] - train_x_i
            kernel_value[i] = np.exp(diff * diff.T / (-2.0 * sigma**2))
    else: # 不使用核函数
        kernel_value = train_x * train_x_i.T
    return kernel_value


def calc_kernel(train_x, kernel_option):
    '''计算核函数矩阵
    input:  train_x(mat):训练样本的特征值
            kernel_option(tuple):核函数的类型以及参数
    output: kernel_matrix(mat):样本的核函数的值
    '''
    m = np.shape(train_x)[0] # 样本的个数
    kernel_matrix = np.mat(np.zeros((m, m))) # 初始化样本之间的核函数值
    for i in range(m):
        kernel_matrix[:, i] = cal_kernel_value(train_x, train_x[i, :], kernel_option)
    return kernel_matrix

def cal_error(svm, alpha_k):
    '''误差值的计算
    input:  svm:SVM模型
            alpha_k(int):选择出的变量
    output: error_k(float):误差值
    '''
    output_k = float(np.multiply(svm.alphas, svm.train_y).T * svm.kernel_mat[:, alpha_k] + svm.b)
    error_k = output_k - float(svm.train_y[alpha_k])
    return error_k


def update_error_tmp(svm, alpha_k):
    '''重新计算误差值
    input:  svm:SVM模型
            alpha_k(int):选择出的变量
    output: 对应误差值
    '''
    error = cal_error(svm, alpha_k)
    svm.error_tmp[alpha_k] = [1, error]

def select_second_sample_j(svm, alpha_i, error_i):
    '''选择第二个样本
    input:  svm:SVM模型
            alpha_i(int):选择出的第一个变量
            error_i(float):E_i
    output: alpha_j(int):选择出的第二个变量
            error_j(float):E_j
    '''
    # 标记为已被优化
    svm.error_tmp[alpha_i] = [1, error_i]
    candidateAlphaList = np.nonzero(svm.error_tmp[:, 0].A)[0]
    
    maxStep = 0
    alpha_j = 0
    error_j = 0

    if len(candidateAlphaList) > 1:
        for alpha_k in candidateAlphaList:
            if alpha_k == alpha_i: 
                continue
            error_k = cal_error(svm, alpha_k)
            if abs(error_k - error_i) > maxStep:
                maxStep = abs(error_k - error_i)
                alpha_j = alpha_k
                error_j = error_k
    else: # 随机选择          
        alpha_j = alpha_i
        while alpha_j == alpha_i:
            alpha_j = int(np.random.uniform(0, svm.n_samples))
        error_j = cal_error(svm, alpha_j)
    
    return alpha_j, error_j

def choose_and_update(svm, alpha_i):
    '''判断和选择两个alpha进行更新
    input:  svm:SVM模型
            alpha_i(int):选择出的第一个变量
    '''
    error_i = cal_error(svm, alpha_i) # 计算第一个样本的E_i
    
    # 判断选择出的第一个变量是否违反了KKT条件
    if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\
        (svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):

        # 1、选择第二个变量
        alpha_j, error_j = select_second_sample_j(svm, alpha_i, error_i)
        alpha_i_old = svm.alphas[alpha_i].copy()
        alpha_j_old = svm.alphas[alpha_j].copy()

        # 2、计算上下界
        if svm.train_y[alpha_i] != svm.train_y[alpha_j]:
            L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])
            H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])
        else:
            L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)
            H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])
        if L == H:
            return 0

        # 3、计算eta
        eta = 2.0 * svm.kernel_mat[alpha_i, alpha_j] - svm.kernel_mat[alpha_i, alpha_i] \
                  - svm.kernel_mat[alpha_j, alpha_j]
        if eta >= 0:
            return 0

        # 4、更新alpha_j
        svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta

        # 5、确定最终的alpha_j
        if svm.alphas[alpha_j] > H:
            svm.alphas[alpha_j] = H
        if svm.alphas[alpha_j] < L:
            svm.alphas[alpha_j] = L

        # 6、判断是否结束      
        if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:
            update_error_tmp(svm, alpha_j)
            return 0

        # 7、更新alpha_i
        svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \
                                * (alpha_j_old - svm.alphas[alpha_j])

        # 8、更新b
        b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
                                                    * svm.kernel_mat[alpha_i, alpha_i] \
                             - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
                                                    * svm.kernel_mat[alpha_i, alpha_j]
        b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
                                                    * svm.kernel_mat[alpha_i, alpha_j] \
                             - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
                                                    * svm.kernel_mat[alpha_j, alpha_j]
        if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):
            svm.b = b1
        elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):
            svm.b = b2
        else:
            svm.b = (b1 + b2) / 2.0

        # 9、更新error
        update_error_tmp(svm, alpha_j)
        update_error_tmp(svm, alpha_i)

        return 1
    else:
        return 0

def SVM_training(train_x, train_y, C, toler, max_iter, kernel_option = ('rbf', 0.431029)):
    '''SVM的训练
    input:  train_x(mat):训练数据的特征
            train_y(mat):训练数据的标签
            C(float):惩罚系数
            toler(float):迭代的终止条件之一
            max_iter(int):最大迭代次数
            kerner_option(tuple):核函数的类型及其参数
    output: svm模型
    '''
    # 1、初始化SVM分类器
    svm = SVM(train_x, train_y, C, toler, kernel_option)
    
    # 2、开始训练
    entireSet = True
    alpha_pairs_changed = 0
    iteration = 0
    
    while (iteration < max_iter) and ((alpha_pairs_changed > 0) or entireSet):
        print ("\t iterration: ", iteration)
        alpha_pairs_changed = 0

        if entireSet:
            # 对所有的样本
            for x in range(svm.n_samples):
                alpha_pairs_changed += choose_and_update(svm, x)
            iteration += 1
        else:
            # 非边界样本
            bound_samples = []
            for i in range(svm.n_samples):
                if svm.alphas[i,0] > 0 and svm.alphas[i,0] < svm.C:
                    bound_samples.append(i)
            for x in bound_samples:
                alpha_pairs_changed += choose_and_update(svm, x)
            iteration += 1
        
        # 在所有样本和非边界样本之间交替
        if entireSet:
            entireSet = False
        elif alpha_pairs_changed == 0:
            entireSet = True

    return svm

def svm_predict(svm, test_sample_x):
    '''利用SVM模型对每一个样本进行预测
    input:  svm:SVM模型
            test_sample_x(mat):样本
    output: predict(float):对样本的预测
    '''
    # 1、计算核函数矩阵
    kernel_value = cal_kernel_value(svm.train_x, test_sample_x, svm.kernel_opt)
    # 2、计算预测值
    predict = kernel_value.T * np.multiply(svm.train_y, svm.alphas) + svm.b
    return predict

def cal_accuracy(svm, test_x, test_y):
    '''计算预测的准确性
    input:  svm:SVM模型
            test_x(mat):测试的特征
            test_y(mat):测试的标签
    output: accuracy(float):预测的准确性
    '''
    n_samples = np.shape(test_x)[0] # 样本的个数
    correct = 0.0
    for i in range(n_samples):
        # 对每一个样本得到预测值
        predict=svm_predict(svm, test_x[i, :])
        # 判断每一个样本的预测值与真实值是否一致
        if np.sign(predict) == np.sign(test_y[i]):
            correct += 1
    accuracy = correct / n_samples
    return accuracy

def save_svm_model(svm_model, model_file):
    '''保存SVM模型
    input:  svm_model:SVM模型
            model_file(string):SVM模型需要保存到的文件
    '''
    with open(model_file, 'wb') as f:
        pickle.dump(svm_model,f)