支持向量机part2

8.3 非线性支持向量机与核函数

非线性数据集则使用非线性支持向量机，其主要特点是利用核技巧(kernel trick)。为此，先要介绍核技巧。

8.3.1 核技巧

1. 非线性分类问题

   非线性分类问题是指通过利用非线性模型才能很好地进行分类的问题。

   如下图所示：

   

   一般来说，对给定的一个训练数据集$T = \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\}$，其中，实例$x_i$属于输入空间，$x_i \in \mathcal{X} = R^n$，对应的标记有两类$y_i \in \mathcal{Y} = \{-1,+1\}$，i=1,2,…,N。如果能用$R^n$中的一个超曲面将正负例正确分开，则称这个问题为非线性可分问题。

   非线性问题的解决思路:对非线性数据集进行非线性变换，将非线性问题变换为线性问题，通过解变换后的线性问题的方法求解原来的非线性问题。

   对图8.7所示的例子，通过变换，将左图中的椭圆变换成右图中的直线，将非线性分类问题变换为线性分类问题：

   设原空间为$\mathcal{X} \subset R^2,x=(x^{(1)},x^{(2)})^{\top} \in \mathcal{X}$，新空间为$\mathcal{Z} \subset R^2,z=(z^{(1)},z^{(2)})^{\top} \in \mathcal{Z}$，定义从原空间到新空间的变换(映射):

   $$z = \phi(x) = ((x^{(1)})^2,(x^{(2)})^2)^{\top}$$

   经过变换$z=\phi(x)$，原空间$\mathcal{X} \subset R^2$变换为新空间$\mathcal{Z} \subset R^2$，原空间中的点相应地变换为新空间中的点，原空间中的椭圆:

   $$w_1(x^{(1)})^2 + w_2(x^{(2)})^2 + b = 0$$

   变换成为新空间中的直线:

   $$w_1 z^{(1)}+w_2 z^{(2)} + b = 0$$

   在变换后的新空间里，直线$w_1z^{(1)}+w_2z^{(2)}+b=0$可以将变换后的正负实例点正确分开。

   这样，原空间的非线性可分问题就变成了新空间的线性可分问题。

2. 核函数的定义

   定义8.6(核函数) 设$\mathcal{X}$是输入空间(欧式空间$R^n$的子集或离散集合)，又设$\mathcal{H}$为特征空间(希尔伯特空间)，如果存在一个从$\mathcal{X}$到$\mathcal{H}$的映射:

   $$\phi(x): \mathcal{X} \rightarrow \mathcal{H} \ (8.65)$$

   使得对所有$x,z \in \mathcal{X}$，函数$K(x,z)$满足条件:

   $$K(x,z) = \phi(x)\bullet \phi(z) \ (8.66)$$

   则称$K(x,z)$为核函数，$\phi(x)$为映射函数。

   核技巧的思路:在学习与预测中只定义核函数$K(x,z)$，而不显式地定义映射函数$\phi$。通常，直接计算$K(x,z)$比较容易，而通过$\phi(x)$和$\phi(z)$计算$K(x,z)$并不容易。

3. 核技巧在支持向量机中的应用

   在线性支持向量机的对偶问题的目标函数(8.37)中的内积$x_i.x_j$可以用核函数$K(x_i,x_j)=\phi(x_i)\bullet \phi(x_j)$来代替。此时对偶问题的目标函数为：

   $$W(\alpha) = \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N \alpha_i \alpha_j y_i y_j K(x_i,x_j)-\sum_{i=1}^N \alpha_i \ (8.67)$$

   同样，分类决策函数中的内积也可以用核函数代替，而分类决策函数式(8.56)变成:

   $$f(x) = \mathrm{sign}\left(\sum_{i=1}^{N_s} \alpha_i^* y_i \phi(x_i)\bullet \phi(x)+b^* \right) = \mathrm{sign} \left(\sum_{i=1}^{N_s} \alpha_i^* y_i K(x_i,x)+b^*\right) \ (8.68)$$

   等价于经过映射函数$\phi$将原来的输入空间变换到一个新的特征空间，将输入空间中的内积$x_i\bullet x_j$变换为特征空间中的内积$\phi(x_i)\bullet \phi(x_j)$。

   在新的特征空间中学习线性支持向量机，当映射函数是非线性函数时，学习到的含有核函数的支持向量机是非线性分类模型。

8.3.2 正定核

函数$K(x,z)$满足什么条件才能成为核函数？(通常所说的核函数就是正定核函数)。

假设$K(x,z)$是定义在$\mathcal{X} \times \mathcal{X}$上的对称函数，并且对任意的$x_1,x_2,...,x_m \in \mathcal{X}$，$K(x,z)$关于$x_1,x_2,...,x_m$的Gram矩阵是半正定的。可以依据函数$K(x,z)$，构成一个希尔伯特空间，其步骤是:

1. 首先定义映射$\phi$并构成向量空间S；
2. 然后在S上定义内积构成内积空间
3. 最后将S完备化构成希尔伯特空间

证明如下所示:

1.定义映射，构成向量空间S

先定义映射

$$\phi: x \rightarrow K(.,x) \ (8.69)$$

根据这一映射，对任意$x_i \in \mathcal{X}$，$\alpha_i \in R$，i=1,2,…,m，定义线性组合

$$f(.) = \sum_{i=1}^m \alpha_i K(.,x_i) \ (8.70)$$

考虑由线性组合为元素的集合$\mathcal{S}$。由于集合$\mathcal{S}$对加法和数乘运算是封闭的，所以$\mathcal{S}$构成一个向量空间。

2.在$\mathcal{S}$上定义内积，使其成为内积空间

在$\mathcal{S}$上定义一个运算*:对任意$f,g \in \mathcal{S}$，

$$f(.) = \sum_{i=1}^m \alpha_i K(.,x_i) \ (8.71)$$ $$g(.) = \sum_{j=1}^l \beta_j K(.,z_j) \ (8.72)$$

定义运算*

$$f*g = \sum_{i=1}^m \sum_{j=1}^l \alpha_i \beta_j K(x_i,z_j) \ (8.73)$$

证明运算*是空间$\mathcal{S}$的内积。为此要证:

$$(cf)*g = c(f*g),c \in R \tag{8.74}$$ $$(f+g)*h = f*h+g*h,h\in\mathcal{S} \tag{8.75}$$ $$f*g = g*f \tag{8.76}$$ $$f*f \ge 0 \tag{8.77}$$ $$f*f = 0 \Leftrightarrow f=0 \ (8.78)$$

其中，(8.74)~(8.76)由式(8.70)~式(8.72)及$K(x,z)$的对称性容易得到。

现证式(8.77)。由式(8.70)及式(8.73)可得:

$$f*f = \sum_{i,j=1}^m \alpha_i\alpha_jK(x_i,x_j)$$

由Gram矩阵的半正定性知上式右端非负，即$f*f \ge 0$。

再证式(8.78)。充分性显然。为证必要性，首先证明不等式:

$$|f*g|^2 \le (f*f)(g*g) \ (8.79)$$

设$f,g \in \mathcal{S},\lambda \in R$，则$f+\lambda g \in \mathcal{S}$，于是，

$$(f+\lambda g)* (f+\lambda g) \ge 0$$ $$f*f + 2\lambda(f*g)+\lambda^2(g*g) \ge 0$$

其左端是$\lambda$的二次三项式，非负，其判别式小于等于0，即

$$(f*g)^2 - (f*f)(g*g) \le 0$$

于是式(8.79)得证。现证若$f*f = 0$，则$f=0$。事实上，若

$$f(.) = \sum_{i=1}^m \alpha_i K(.,x_i)$$

则按运算*的定义式(8.73)，对任意的$x \in \mathcal{X}$，有

$$K(.,x)*f = \sum_{i=1}^m \alpha_i K(x,x_i) = f(x)$$

于是，

$$|f(x)|^2 = |K(.,x)*f|^2 \ (8.80)$$

由式(8.79)和式(8.77)有:

$$|K(.,x)*f|^2 \le (K(.,x)*K(.,x))(f*f) = K(x,x)(f*f)$$

由式(8.80)有

$$|f(x)|^2 \le K(x,x)(f*f)$$

此式表明，当$f*f = 0$时，对任意的x都有$|f(x)|=0$。

至此，证明了为向量空间$\mathcal{S}$的内积。赋予内积的向量空间为内积空间。*因此$\mathcal{S}$是一个内积空间。

既然*为$\mathcal{S}$的内积运算，那么仍然用$\bullet$表示，即若:

$$f(.) = \sum_{i=1}^m \alpha_i K(.,x_i)\ , \ g(.) = \sum_{i=1}^l \beta_j K(.,z_j)$$

则，

$$f\bullet g = \sum_{i=1}^m \sum_{j=1}^l \alpha_i \beta_j K(x_i,z_j) \ (8.81)$$

3.将内积空间$\mathcal{S}$完备化为希尔伯特空间

由式(8.81)定义的内积可以得到范数:

$$\|f\| = \sqrt{f \bullet f} \ (8.82)$$

因此，$\mathcal{S}$是一个赋范向量空间。根据泛函分析理论，对于不完备的赋范向量空间$\mathcal{S}$，一定可以使之完备化，得到完备的赋范向量空间$\mathcal{H}$。

一个内积空间，当作为一个赋范向量空间是完备的时候，就是希尔伯特空间。

核函数K具有再生性，即满足:

$$K(.,x)\bullet f = f(x) \ (8.83)$$

及

$$K(.,x)\bullet K(.,z) = K(x,z) \ (8.84)$$

称为再生核。

4.正定核的充要条件

定理8.5(正定核的充要条件) 设$K: \mathcal{X} \times \mathcal{X} \rightarrow R$是对称函数，则$K(x,z)$为正定核函数的充要条件是对任意$x_i \in \mathcal{X}$，i=1,2,…,m, $K(x,z)$对应的Gram矩阵：

$$K = \left[K(x_i,x_j)\right]_{m\times m} \ (8.85)$$

是半正定矩阵。

证明:

必要性。由于$K(x,z)$是$\mathcal{X} \times \mathcal{X}$上的正定核，所以存在从$\mathcal{X}$到希尔伯特空间$\mathcal{H}$的映射$\phi$，使得:

$$K(x,z) = \phi(x) \bullet \phi(z)$$

于是，对任意$x_1,x_2,...,x_m$，构造$K(x,z)$关于$x_1,x_2,...,x_m$的Gram矩阵:

$$\left[K_{ij}\right]_{m\times m} = \left[K(x_i,x_j)\right]_{m \times m}$$

对任意$c_1,c_2,...,c_m \in R$，有

$$\begin{aligned}\sum^{m}_{i,j=1} c_ic_j K(x_i,x_j) &= \sum_{i,j=1}^m c_ic_j(\phi(x_i)\bullet \phi(x_j)) \\ &= \left(\sum_i c_i \phi(x_i)\right) \bullet \left(\sum_j c_j \phi(x_j)\right) = \left\|\sum_i c_i \phi(x_i)\right\|^2 \ge 0\end{aligned}$$

表明$K(x,z)$关于$x_1,x_2,...,x_m$的Gram矩阵是半正定的。

充分性。已知对称函数$K(x,z)$对任意$x_1,x_2,...,x_m \in \mathcal{X}$，$K(x,z)$关于$x_1,x_2,...,x_m$的Gram矩阵是半正定的(即$x_i K(x,z) x_i \ge 0$)。根据前面的结果，对给定的$K(x,z)$，可以构造从$\mathcal{X}$到某个希尔伯特空间$\mathcal{H}$的映射:

$$\phi:x \rightarrow K(.,x) \ (8.86)$$

由式(8.83)可知，

$$K(.,x)\bullet f = f(x)$$

并且

$$K(.,x)\bullet K(.,z) = K(x,z)$$

由式(8.86)即得

$$K(x,z) = \phi(x) \bullet \phi(z)$$

表明$K(x,z)$是$\mathcal{X} \times \mathcal{X}$上的核函数。

定义8.7（正定核的等价定义） 设$\mathcal{X} \subset R^n$，$K(x,z)$是定义在$\mathcal{X} \times \mathcal{X}$上的对称函数，如果对任意$x_i \in \mathcal{X}$，i=1,2,…,m，$K(x,z)$对应的Gram矩阵

$$K = \left[K(x_i,x_j)\right]_{m \times m} \ (8.87)$$

是半正定矩阵，则称$K(x,z)$是正定核。

正定核的检验比较难，因为要求对任意$x_i \in \mathcal{X}$，i=1,2,…,m，验证K对应的Gram矩阵是否为半正定的。

所以往往使用已有的核函数。另外，由Mercer定义可以得到Mercer核，正定核比Mercer核更具一般性。

8.3.3 常用核函数

1.多项式核函数

$$K(x,z) = (x.z+1)^p \ (8.88)$$

对应的支持向量机是一个p次多项式分类器。在此情形下，分类决策函数成为

$$f(x) = \mathrm{sign} \left(\sum^{N_s}_{i=1} \alpha_i^* y_i (x_i \bullet x+1)^p+b^*\right) \ (8.89)$$

2.高斯核函数

$$K(x,z) = \exp \left(-\frac{\|x-z\|^2}{2\sigma^2}\right)\ (8.90)$$

对应的支持向量机是高斯径向基函数分类器。在此情形下，分类决策函数成为:

$$f(x) = \mathrm{sign} \left( \sum_{i=1}^{N_s} \alpha_i^* y_i \exp \left(- \frac{\|x-z\|^2}{2\sigma^2}\right) + b^*\right) \ (8.91)$$

3.字符串核函数

核函数不仅仅可以定义在欧式空间上，还可以定义在离散数据的集合上。比如，字符串核是定义在字符串集合上的核函数。字符串核函数在文本分类、信息检索、生物信息学等方面都有应用。

字符串核函数$k_n(s,t)$给出了字符串s和t中长度等于n的所有子串组成的特征向量的余弦相似度。直观上，两个字符串相同的子串越多，它们就越相似，字符串核函数的值就越大。字符串核函数可以由动态规划快速地计算。

8.3.4 非线性支持向量分类机

如上所述，利用核技巧，可以将线性分类的学习方法应用到非线性分类问题中去。将线性支持向量机扩展到非线性支持向量机，只需将线性支持向量机对偶形式中的内积换成核函数。

定义8.8(非线性支持向量机) 从非线性分类训练集，通过核函数与软间隔最大化，或凸二次规划(8.95~8.97)，学习得到的分类决策函数。

$$f(x) = \mathrm{sign} \left(\sum_{i=1}^N \alpha_i^* y_i K(x,x_i)+b^*\right) \ (8.94)$$

称为非线性支持向量，$K(x,z)$是正定核函数。

算法8.4（非线性支持向量机学习算法）

输入:训练数据集$T = \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\}$，其中$x_i \in \mathcal{X} = R^n$，$y_i \in \mathcal{Y} = \{-1,+1\}$，i=1,2,…,N;

输出:分类决策函数

(1)选取适当的核函数$K(x,z)$和适当的参数C，构造并求解最优化问题:

$$\min_{\alpha} \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N \alpha_i \alpha_j y_i y_j K(x_i,x_j)-\sum_{i=1}^N \alpha_i \tag{8.95}$$ $$\mathrm{s.t.} \sum_{i=1}^N \alpha_i y_i = 0 \tag{8.96}$$ $$0 \le \alpha_i \le C，i=1,2,…,N \tag{8.97}$$

求得最优解$\alpha^* = (\alpha_1^*,\alpha_2^*,...,\alpha_N^*)^{\top}$。

(2)选择$\alpha^*$的一个正分量$0 < \alpha_j^* < C$，计算

$$b^* = y_j - \sum_{i=1}^N \alpha_i^* y_i K(x_i \bullet x_j)$$

(3)构造决策函数:

$$f(x) = \mathrm{sign} \left(\sum_{i=1}^N \alpha_i^* y_i K(x\bullet x_i)+b^*\right)$$

当$K(x,z)$是正定核函数时，问题(8.95)~(8.97)是凸二次规划问题，解是存在的。

8.4 序列最小最优化算法

SMO(sequential minimal optimization)算法要解如下凸二次规划的对偶问题:

$$\min_{\alpha} \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N \alpha_i \alpha_j y_i y_j K(x_i,x_j)-\sum_{i=1}^N \alpha_i \tag{8.98}$$ $$\mathrm{s.t.} \sum_{i=1}^N \alpha_iy_i = 0 \tag{8.99}$$ $$0\le \alpha_i \le C，i=1,2,...,N \tag{8.100}$$

在这个问题中，变量是拉格朗日乘子，一个变量$\alpha_i$对应于一个样本点$(x_i,y_i)$；变量的总数等于训练样本容量N。

SMO算法包括两个部分：求解两个变量二次规划的解析方法和选择变量的启发式方法。

8.4.1两个变量二次规划的求解方法

不失一般性，假设选择的两个变量是$\alpha_1$，$\alpha_2$，其他变量$\alpha_i(i=3,4,...,N)$是固定的。于是SMO的最优化问题(8.98)~(8.100)的子问题可以写成:

$$\begin{aligned}\min_{\alpha_1,\alpha_2} W(\alpha_1,\alpha_2) = &\frac{1}{2} K_{11}\alpha_1^2 + \frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2 \\ &-(\alpha_1+\alpha_2)+y_1\alpha_1 \sum_{i=3}^N y_i \alpha_i K_{il}+y_2\alpha_2\sum_{i=3}^N y_i \alpha_i K_{i2}\end{aligned} \ (8.101)$$ $$\mathrm{s.t.} \alpha_1y_1 + \alpha_2y_2 = -\sum_{i=3}^N y_i \alpha_i = \zeta \ (8.102)$$ $$0\le \alpha_i \le C \ , \ i=1,2 \ (8.103)$$

其中，$K_{ij} = K(x_i,x_j)$，i,j=1,2,…,N，$\zeta$是常数，目标函数式(8.101)中省略了不含$\alpha_1$，$\alpha_2$的常数项。

为了求解两个变量的二次规划问题(8.101)~(8.103)：

1. 首先分析约束条件
2. 然后在此约束条件下求极小。

8.4.2 分析约束条件

由于只有两个变量$(\alpha_1,\alpha_2)$，约束可以用二维空间中的图形表示，如下图所示:

不等式约束(8.103)使得$(\alpha_1,\alpha_2)$在盒子$[0,C]\times [0,C]$内，等式约束(8.102)使$(\alpha_1,\alpha_2)$在平行于盒子$[0,C]\times [0,C]$的对角线的直线上。因此要求的是目标函数在一条平行于对角线的线段上的最优值。

假设问题(8.101)~(8.103)的初始可行解为$\alpha_1^{old}$，$\alpha_2^{old}$，最优解为$\alpha_1^{new}$，$\alpha_2^{new}$，并且假设在沿着约束方向未经剪辑时$\alpha_2$的最优解为$\alpha_2^{new,unc}$。

由于$\alpha_2^{new}$满足不等式约束条件(8.103)，所以最优值$\alpha_2^{new}$的取值范围必须满足条件:

$$L \le \alpha_2^{new} \le H$$

其中，L与H是$\alpha_2^{new}$所在的对角线段端点的界。

如果$y_1 \ne y_2$，则

$$L = \max(0,\alpha_2^{old}-\alpha_1^{old})$$ $$H = \min(C,C+\alpha_2^{old}-\alpha_1^{old})$$

L式子max函数的第二项，对应于$\alpha_1^{new}=0$，$\alpha_2^{new}=-k$，而$-k=\alpha_2^{old}-\alpha_1^{old}$

H式子max函数的第二项，对应于$\alpha_1^{new}=C$，同理可得以上表达式。

如果$y_1=y_2$，则

$$L = \max(0, \alpha_2^{old}+\alpha_1^{old}-C)$$ $$H = \min(C, \alpha_2^{old}+\alpha_1^{old})$$

首先，求沿着约束方向未经剪辑即未考虑不等式约束(8.103)时$\alpha_2$的最优解$\alpha_2^{new,unc}$；

然后，再求剪辑后$\alpha_2$的解$\alpha_2^{new}$。我们用定理来叙述这个结果，为了叙述简单，记:

$$g(x) = \sum_{i=1}^N \alpha_i y_i K(x_i,x)+b \tag{8.104}$$

令

$$E_i = g(x_i)-y_i = \left(\sum_{j=1}^N \alpha_j y_j K(x_j,x_i)+b\right) - y_i \ , i=1,2 \tag{8.105}$$

当i=1,2时，$E_i$为函数g(x)对输入$x_i$的预测值与真实输出$y_i$之差。

定理8.6 最优化问题(8.101)~(8.103)沿着约束方向未经剪辑时的解是:

$$\alpha_2^{new,unc} = \alpha_2^{old} + \frac{y_2(E_1-E_2)}{\eta} \tag{8.106}$$

其中，

$$\eta = K_{11}+K_{22}-2K_{12} = \|\Phi(x_1) -\Phi(x_2)\|^2 \tag{8.107}$$

关于$\Phi(x)$是输入空间到特征空间的映射，$E_i$，i=1,2，由式(8.105)给出。

经剪辑后$\alpha_2$的解是

$$\alpha_2^{new} = \left\{\begin{aligned}&H, & \alpha_2^{new,unc} > H \\ &\alpha_2^{new,unc},&L\le \alpha_2^{new,unc} \le H \\ &L,&\alpha_2^{new,unc}<L\end{aligned} \right . \ (8.108)$$

由$\alpha_2^{new}$求得$\alpha_1^{new}$是(根据式(8.102))

$$\alpha_1^{new} = \alpha_1^{old} + y_1y_2(\alpha_2^{old}-\alpha_2^{new}) \ (8.109)$$

证明:

引进记号：

$$v_i = \sum_{j=3}^{N}\alpha_jy_jK_{ij}=g(x_i)-\sum_{j=1}^{2}\alpha_jy_jK_{ij}-b \ , i=1,2$$

目标函数可写成

$$W(\alpha_1,\alpha_2) = \frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2-(\alpha_1+\alpha_2)+y_1v_1\alpha_1+y_2v_2\alpha_2 \ (8.110)$$

由$\alpha_1y_1 = \zeta -\alpha_2y_2$及$y_i^2=1$，可将$\alpha_1$表示为：

$$\alpha_1 = (\zeta-\alpha_2y_2)y_1$$

代入式(8.110)，得到只是$\alpha_2$的函数的目标函数:

$$W(\alpha_2)=\frac{1}{2}K_{11}( \zeta-\alpha_2y_2)^2+\frac{1}{2}K_{22}\alpha_2^2+y_2K_{12}(\zeta-\alpha_2y_2)\alpha_2- (\zeta-\alpha_2y_2)y_1-\alpha_2+v_1 (\zeta-\alpha_2y_2)+y_2v_2\alpha_2$$

对$\alpha_2$求导数

$$\frac{\partial W}{\partial \alpha_2} = K_{11}\alpha_2 + K_{22}\alpha_2-2K_{12}\alpha_2-K_{11}\zeta y_2+K_{12}\zeta y_2+y_1y_2-1-v_1y_2+y_2v_2$$

令其为0，得到

$$\begin{aligned}(K_{11}+K_{22}-2K_{12})\alpha_2 = &y_2(y_2-y_1+\zeta K_{11}-\zeta K_{12}+v_1-v_2) \\ = &y_2 \left[y_2-y_1+\zeta K_{11}-\zeta K_{12}+\left(g(x_1)-\sum_{j=1}^2 y_j\alpha_j K_{1j}-b\right)\right] \\ &-\left(g(x_2)-\sum_{j=1}^2 y_2 \alpha_j K_{2j}-b\right)\end{aligned}$$

将$\zeta = \alpha_1^{old}y_1+\alpha_2^{old}y_2$代入，得到:

$$\begin{aligned}(K_{11}+K_{22}-2K_{12})\alpha_2^{new,unc} &= y_2((K_{11}+K_{22}-2K_{12})\alpha_2^{old}y_2+y_2-y_1+g(x_1)-g(x_2)) \\ &= (K_{11}+K_{22}-2K_{12})\alpha_2^{old}+y_2(E_1-E_2) \end{aligned}$$

将$\eta = K_{11}+K_{22}-2K_{12}$代入，于是得到：

$$\alpha_2^{new,unc} = \alpha_2^{old} + \frac{y_2(E_1-E_2)}{\eta}$$

要使其满足不等式约束必须将其限制在区间[L,H]内，从而得到$\alpha_2^{new}$的表达式(8.108)。由等式约束(8.102)得到$\alpha_1^{new}$的表达式(8.109)。于是得到最优化问题(8.101)~(8.103)的解($\alpha_1^{new}$，$\alpha_2^{new}$)。

8.4.3 变量的选择方法

SMO算法在每个子问题中选择两个变量优化，其中至少一个变量是违反KKT条件的。

1.第 1个变量的选择

SMO称选择第1个变量的过程为外层循环。外层循环在训练样本中选取违反KKT条件最严重的样本点，并将其对应的变量作为第1个变量。具体地，检验训练样本点$(x_i,y_i)$是否满足KKT条件，即

$$\alpha_i = 0 \Leftrightarrow y_ig(x_i) \ge 1 \ (8.111)$$ $$0 < \alpha_i < C \Leftrightarrow y_ig(x_i) = 1 \ (8.112)$$ $$\alpha_i = C \Leftrightarrow y_ig(x_i) \le 1 \ (8.113)$$

其中，$g(x_i)=\sum_{j=1}^N \alpha_jy_jK(x_i,x_j)+b$。

该检验是在$\varepsilon$指定阈值范围内进行的。在检验过程中，外层循环首先遍历所有满足条件$0<\alpha_i<C$的样本点，即在间隔边界上的支持向量点，检验它们是否满足KKT条件。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验它们是否满足KKT条件。

2.第2个变量的选择

SMO称选择第2个变量的过程为内层循环。假设在外层循环中已经找到第1个变量$\alpha_1$，现在要在内层循环中找第2个变量$\alpha_2$。第2个变量选择的标准是希望能使$\alpha_2$有足够大的变化。

由式(8.106)和式(8.108)可知，$\alpha_2^{new}$是依赖于$|E_1-E_2|$的。一种简单做法是选择$\alpha_2$，使其对应的$|E_1-E_2|$最大。(因为$\alpha_1$已定，$E_1$也确定了。如果$E_1$是正的，那么选择最小的$E_i$作为$E_2$;如果$E_1$是负的，那么选择最大的$E_i$作为$E_2$)。为节省时间，把所有$E_i$保存到一个列表中。

在特殊情况下，如果内层循环通过以上方法选择的$\alpha_2$不能使目标函数有足够的下降，那么采用以下启发式规则继续选择$\alpha_2$。

启发式规则如下所示：

1. 遍历在间隔边界上的支持向量点，依次将其对应的变量作为$\alpha_2$使用，直到目标函数有足够的下降(使用指定阈值$\varepsilon$来标识)。
2. 若找不到合适的$\alpha_2$，那么遍历训练数据集
3. 若仍找不到合适的$\alpha_2$，则放弃第1个$\alpha_1$，再通过外层循环寻求另外的$\alpha_1$。

3.计算阈值b和差值$E_i$

在每次完成两个变量的优化后，都要重新计算阈值b。当$0 < \alpha_1^{new} <C$时，由KKT条件(8.112)可知：

$$\sum_{i=1}^N \alpha_i y_i K_{i1}+b = y_1$$

于是，

$$b_1^{new} = y_1 - \sum_{i=3}^N \alpha_i y_i K_{i1}-\alpha_1^{new}y_1K_{11}-\alpha_2^{new}y_2 K_{21} \tag{8.114}$$

由$E_1$的定义式(8.105)有

$$E_1 = \sum_{i=3}^N \alpha_i y_i K_{i1} + \alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old}-y_1$$

式(8.114)的前两项可写成：

$$y_1 - \sum_{i=3}^N \alpha_i y_i K_{i1} = -E_1+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old}$$

代入式(8.114)，可得：

$$b_1^{new} = -E_1-y_1K_{11}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{21}(\alpha_2^{new}-\alpha_2^{old})+b^{old} \tag{8.115}$$

同样，如果$0 < \alpha_2^{new}<C$，那么，

$$b_2^{new} = -E_2-y_1K_{12}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{22}(\alpha_2^{new}-\alpha_2^{old})+b^{old} \tag{8.116}$$

如果$\alpha_1^{new}$，$\alpha_2^{new}$同时满足条件$0<\alpha_i^{new}<C$，i=1,2，那么$b_1^{new}=b_2^{new}$。

如果$\alpha_1^{new}$，$\alpha_2^{new}$是0或者C，那么$b_1^{new}$和$b_2^{new}$以及它们之间的数都是符合KKT条件的阈值，这时选择他们的中点作为$b^{new}$。

在每次完成两个变量的优化之后，还必须更新对应的$E_i$值，并将它们保存在列表中。

关于$E_i$值的更新要用到$b^{new}$值，以及所有支持向量对应的$\alpha_j$：

$$E_i^{new} = \sum_{S} y_j \alpha_j K(x_i,x_j)+b^{new} - y_i \tag{8.117}$$

其中,S是所有支持向量$x_j$的集合。

8.4.4 SMO算法

算法8.5(SMO算法)

输入：训练数据集$T = \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\}$，其中，$x_i \in \mathcal{X}=R^n$，$y_i \in \mathcal{Y}=\{-1,+1\}$，i=1,2,…,N，精度$\varepsilon$；

输出：近似解$\hat{\alpha}$

1. 取初值$\alpha^{(0)}=0$，令k=0

2. 选取优化变量$\alpha_1^{(k)}$，$\alpha_2^{(k)}$，解析求解两个变量的最优化问题(8.101)~(8.103)，求得最优解$\alpha_1^{(k+1)}$，$\alpha_2^{(k+1)}$，更新$\alpha$为$\alpha^{(k+1)}$；

3. 若在精度$\varepsilon$范围内满足停机条件

   $$\sum_{i=1}^N \alpha_i y_i = 0$$ $$0 \le \alpha_i \le C \ , i=1,2,...,N$$ $$y_i \bullet g(x_i) = \left\{ \begin{aligned}&\ge1 , &\{x_i|\alpha_i=0\}\\&=1, &\{x_i|0<\alpha_i<C\}\\ &\le 1, &\{x_i|\alpha_i=C\}\end{aligned} \right.$$

​ 其中，

$$g(x_i) \sum_{j=1}^N \alpha_j y_j K(x_j,x_i)+b$$

​ 则转步骤4；否则令k=k+1,转步骤2

1. 取$\hat{\alpha}=\alpha^{k+1}$。

class MySVM(object):
    def __init__(self, kernel='linear', epsilon=0.00001, C=1000, n_iter=5000):
        self.kernel = kernel #核函数
        self.epsilon = epsilon #用于检验kkt条件的误差范围
        self.C = C #式(8.103)
        self.Max_Iteration=n_iter
        
    def _init_parameters(self, features, labels):
        self.X = features
        self.Y = labels
        
        self.b = 0.0
        self.N,self.n = features.shape
        self.alpha = np.zeros(self.N) #对偶算法学习参数
        self.E = [self._E_(i) for i in xrange(self.N)] #误差
        
    def _satisfy_KKT(self, i):
        ygx = self.Y[i] * self._g_(i)
        if abs(self.alpha[i])<self.epsilon: #公式(8.111)
            return ygx>=1
        elif abs(self.alpha[i]-self.C)<self.epsilon:# 公式(8.113)
            return ygx <= 1
        else: # 公式(7.112)
            return abs(ygx-1)<self.epsilon
    
        
    def is_stop(self):
        """
        所有的样例都满足KKT条件
        """
        result = 0
        for in xrange(self.N):
            result += self.alpha[i]*self.Y[i]
        if result>self.epsilon:
            return False
        
        for i in xrange(self.N):
            satisfy = self._satisfy_KKT(i)
            
            if not satisfy:
                return False
        return True
    
    
    def _select_two_parameters(self):
        #第一部分是alphai满足条件(0,C)的样本点，即间隔边界上的支持向量点
        i1_list_1 = [i for i in xrange(self.N) if self.alpha[i]>0 and self.alpha[i]<self.C]
        #第二部分是当所有支持向量点满足KKT条件时，遍历其他所有非支持向量点
        i1_list_2 = [i for i in xrange(self.N) if self.alpha[i]<=0 or self.alpha[i]>=self.C]
        
        i1_list = i1_list_1
        i1_list.extend(i1_list_2)
        
        for i in i1_list:
            if self._satisfy_KKT(i):
                continue
            #找到一个违反KKT条件的样本点
            E1 = self.E[i]
            max_ = (0,0)
            
            for j in range(self.N): #选择一个使得|Ei-Ej|最大的j
                if j==i:
                    continue
                    
                E2 = self.E[j]
                if abs(E1-E2)>max_[0]:
                    max_ = (abs(E1-E2),j)
                    
            return i,max_[1]
        
    def _K_(self, x1, x2):
        '''
        核函数,这里只实现线性核函数，相当于线性支持向量机
        
        '''
        if self.kernel=='linear':
            return np.dot(x1,x2)
        
        return np.dot(x1,x2)
    
    def _g_(self, i):
        '''公式(7.104)'''
        result = self.b
        for j in xrange(self.N):
            result += self.alpha[j]*self.Y[j]*self._K_(self.X[i],self.X[j])
            
        return result
    
    def _E_(self, i):
        '''公式(7.105)'''
        return self._g_(i) - float(self.Y[i])
    
    def train(self, features, labels):
        self._init_parameters(features, labels)
        
        for i in range(self.Max_Iteration):
            i1, i2 = self._select_two_parameters()
            
            if self.Y[i1]==self.Y[i2]:
                L = max(0, self.alpha[i2]+self.alpha[i1]-self.C)
                H = min(self.C, self.alpha[i2]+self.alpha[i1])
            else:
                #Y[i1]!=self.Y[i2]
                L = max(0, self.alpha[i2]-self.alpha[i1])
                H = min(self.C, self.C+self.alpha[i2]-self.alpha[i1])
                
            E1 = self.E[i1]
            E2 = self.E[i2]
            #公式(7.107)
            eta = self._K_(self.X[i1], self.X[i1])+self._K_(self.X[i2],self.X[i2])\
                  -2*self._K_(self.X[i1],self.X[i2])
            #print eta
            #未经剪辑的alpha2
            alpha2_new_unc = self.alpha[i2]+self.Y[i2]*(E1-E2)/eta #公式(8.106)
            
            #剪辑后的alpha2#公式(8.108)
            if alpha2_new_unc>H:
                alpha2_new = H
            elif alpha2_new_unc<L:
                alpha2_new = L
            else:
                alpha2_new = alpha2_new_unc
            
            #剪辑后的alpha1；公式(8.109)
            alpha1_new = self.alpha[i1]+self.Y[i1]*self.Y[i2]*(self.alpha[i2]-alpha2_new)
            
            #公式(8.115)及公式(8.116)
            b_new = 0
            b1_new = -E1-self.Y[i1]*self._K_(self.X[i1], self.X[i1])*(alpha1_new-self.alpha[i1]) \
                        -self.Y[i2]*self._K_(self.X[i2], self.X[i1])*(alpha2_new-self.alpha[i2])+self.b
            
            b2_new = -E2-self.Y[i1]*self._K_(self.X[i1], self.X[i2])*(alpha1_new-self.alpha[i1]) \
                        -self.Y[i2]*self._K_(self.X[i2], self.X[i2])*(alpha2_new-self.alpha[i2])+self.b
            
            if alpha1_new>0 and alpha1_new<self.C:
                b_new = b1_new
            elif alpha2_new>0 and alpha2_new<self.C:
                b_new = b2_new
            else:
                b_new = (b1_new+b2_new)/2
            
            #更新alpha和b
            self.alpha[i1] = alpha1_new
            self.alpha[i2] = alpha2_new
            self.b = b_new
            #更新Ei
            self.E[i1] = self._E_(i1)
            self.E[i2] = self._E_(i2)
            
        self.w = self.calcWs(features, labels)
        
    def _predict_(self, feature):
        result = self.b
        
        for i in range(self.N):
            result += self.alpha[i]*self.Y[i]*self._K_(feature, self.X[i])
            
        if result>0:
            return 1
        
        return -1
    
    def predict(self, features):
        results = []
        
        for feature in features:
            results.append(self._predict_(feature))
            
        return results
    
    def calcWs(self, X, y):
        '''
        基于alpha计算w值
        '''
        m,n = X.shape
        w = np.zeros(n)

        for i in range(m):
            #print(X[i,:].shape)
            w += X[i,:].T*self.alpha[i]*y[i]
        return w
    
    def plotSVM(self, features, labels):
        b = self.b
        fig = plt.figure()
        ax = fig.add_subplot(111)
        
        uni_label = np.unique(labels)
        
        feature_0 = np.array(features[labels==uni_label[0]])
        feature_1 = np.array(features[labels==uni_label[1]])

        
        ax.scatter(feature_0[:,0],feature_0[:,1],color='blue', marker='o',label=str(uni_label[0]))
        ax.scatter(feature_1[:,0],feature_1[:,1],color='blue', marker='x',label=str(uni_label[1]))
        #决策边界
        x_min = min(features[:,0])
        x_max = max(features[:,0])
        x = np.arange(x_min, x_max, 1.2)
        
        y_min = min(features[:,1])
        y_max = max(features[:,1])
        y = (-b-self.w[0]*x)/self.w[1]
        
        x = x[y<y_max]
        y = y[y<y_max]      
        
        ax.plot(x, y)
        
        #找到支持向量，并在图中标红
        for i in range(len(labels)):
            if self.alpha[i] > 0.0 and labels[i]==uni_label[0]:
                ax.plot(features[i,0], features[i,1], 'ro')
            elif self.alpha[i] > 0.0:
                ax.plot(features[i,0], features[i,1], 'rx')
        plt.legend(loc='upper left')
        plt.show()

#另一个版本的与统计学习方法上有所区别

class SVM(object):
    def __init__(self, C=1000, toler=0.001, n_iter=40):
        self.C = C
        self.toler = toler
        self.n_iter = n_iter
        
    def selectJrand(self,i,m):
        '''
        《统计学习方法》中对j的选择是使|Ei-Ej|最大；但是这里是随机选择一个
        '''
        j = i
        while i==j:
            j = int(np.random.uniform(0,m))
        return j
    
    def clipAlpha(self, aj, H, L):
        '''
        更新aj,公式(7.108)
        '''
        if aj>H:
            return H
        if aj<L:
            return L
        return aj
    
    def _satisfy_KKT(self,i, Yi, gx):
        ygx = Yi * gx
        if abs(self.alphas[i])<self.toler:
            return ygx>=1
        elif abs(self.alphas[i]-self.C)<self.toler:
            return ygx <= 1
        else:
            return abs(ygx-1)<self.toler
    
    def _E_(self, gx, yi):
        #式(7.105)
        return gx-yi
    
    def _K_(self, xi, xj):
        '''
        核函数
        '''
        return np.dot(xi, xj)

    def fit(self, X, y):
        m,n = X.shape
        #初始化b和alphas(alpha有点类似权重值)
        self.b = 0
        self.alphas = np.zeros(m)

        for cur_iter in range(self.n_iter):
            alphaPairsChanged = 0
            
            #遍历所有的样本点，找到第一个不满足KKT条件的样本点作为第1个变量
            for i in range(m):
                #求最优间隔分类器的净输入
                output_i = self.predict_proba(X,y,i)#(式7.104)
                error_i = self._E_(output_i,float(y[i]))
                #判断是否满足KKT条件
                if self._satisfy_KKT(i, y[i], output_i):
                    continue
                
                j = self.selectJrand(i,m)
                output_j = self.predict_proba(X,y,j)
                error_j = self._E_(output_j,float(y[j]))
                
                alphaIold = self.alphas[i]
                alphaJold = self.alphas[j]
                
                if y[i]!=y[j]:
                    L = max(0, self.alphas[j]-self.alphas[i])
                    H = min(self.C, self.C+self.alphas[j]-self.alphas[i])
                else:
                    L = max(0, self.alphas[j]+self.alphas[i]-self.C)
                    H = min(self.C, self.alphas[j]+self.alphas[i])

                    
                #公式(7.107)
                eta = self._K_(X[i,:],X[i,:])+self._K_(X[j,:],X[j,:])-2.0*self._K_(X[i,:],X[j,:])
                #print eta
                if eta < 0:
                    continue
                    
                #未经剪辑的alphaj:更新alphas[j]
                alphaj_new = self.alphas[j]+y[j]*(error_i-error_j)/eta
                
                #经过剪辑的alphaj
                alphaj_new = self.clipAlpha(alphaj_new, H, L)
                self.alphas[j] = alphaj_new
                
                if(abs(self.alphas[j]-alphaJold)<0.00001):
                    continue
                    
                #更新alphas[i]
                self.alphas[i] = alphaIold + y[j]*y[i]*(alphaJold-self.alphas[j])
                
                #更新b
                b1_new = self.b -error_i-y[i]*(self.alphas[i]-alphaIold)*self._K_(X[i,:],X[i,:])\
                                        -y[j]*(self.alphas[j]-alphaJold)*self._K_(X[i,:],X[j,:])
                
                b2_new = self.b-error_j-y[i]*(self.alphas[i]-alphaIold)*self._K_(X[i,:],X[j,:])\
                                       -y[j]*(self.alphas[j]-alphaJold)*self._K_(X[j,:],X[j,:])
                
                b_new = (b1_new+b2_new)/2
                
                self.b = b_new
                alphaPairsChanged += 1
            
            if alphaPairsChanged == 0:
                continue
            else:
                break
                
        self.w = self.calcWs(X,y)
        w_norm = np.linalg.norm(self.w)
        self.w = self.w/w_norm
        self.b = self.b/w_norm
        return self
    
    def predict_proba(self, X, y, i):
        '''
        公式(7.104)
        '''
        result = self.b
        m,n = X.shape
        for j in range(m):
            result += self.alphas[j]*y[j]*self._K_(X[i,:],X[j,:])
        return result
    
    def calcWs(self, X, y):
        '''
        基于alpha计算w值
        '''
        m,n = X.shape
        w = np.zeros(n)

        for i in range(m):
            #print(X[i,:].shape)
            w += X[i,:].T*self.alphas[i]*y[i]
        return w
    
    def plotSVM(self, features, labels):
        features = np.matrix(features)
        #b = np.array(self.b)[0]
        b = self.b
        fig = plt.figure()
        ax = fig.add_subplot(111)
        
        uni_label = np.unique(labels)
        features_1 = features[np.array(labels==uni_label[0]).ravel()]
        features_0 = features[np.array(labels==uni_label[1]).ravel()]
        
        ax.scatter(features_1[:,0].flatten().A[0], features_1[:,1].flatten().A[0], color='blue', marker='o', label=str(uni_label[0]))
        ax.scatter(features_0[:,0].flatten().A[0], features_0[:,1].flatten().A[0], color='blue', marker='x', label=str(uni_label[1]))
        
        #决策边界
        x_min = min(features[:,0])
        x_max = max(features[:,0])
        print(x_min)
        print(x_max)
        
        x = np.arange(x_min, x_max, 0.1)
        y = (-b-self.w[0]*x)/self.w[1]
        ax.plot(x, y)
        
        #找到支持向量，并在图中标红
        for i in range(len(labels)):
            if self.alphas[i] > 0.0 and labels[i]==uni_label[0]:
                ax.plot(features[i,0], features[i,1], 'ro')
            elif self.alphas[i] > 0.0:
                ax.plot(features[i,0], features[i,1], 'rx')
        plt.legend(loc='upper left')
        plt.show()