文章目录

一、原始目标（一）超平面方程（二）严格线性可分问题二、目标解析三、兼容软间隔四、兼容非线性分割五、拉格朗日乘数法六、核函数七、SMO算法（一）无约束最小化（二）根据约束条件做修剪八、应用示例（一）SVM函数代码（二）绘图代码（三）线性分割效果（四）非线性分割效果

通俗来讲，所谓支持向量机是一种分类器，对于做出标记的两组向量，给出一个最优分割超曲面把这两组向量分割到两边，使得两组向量中离此超平面最近的向量（即所谓支持向量）到此超平面的距离都尽可能远。

一、原始目标

（一）超平面方程

所谓向量，可以简单理解成在直角坐标系中由原点出发，指向某一坐标点的一个剪头，因此它只有两个特征，大小和方向，如下图（以二维空间为例）中x\boldsymbol xx。所谓nnn维超平面（hyperplane）是指在nnn维内积空间中，所有在同一个法向量（如下图中w\boldsymbol ww）上投影长度相同的向量组成的集合，如下图中超平面。

设法向量w=[w1w2⋯wn]T\boldsymbol w = [\begin{matrix} w_1 & w_2 & \cdots & w_n\end{matrix}]^Tw=[w1​​w2​​⋯​wn​​]T，其模长为：

∥w∥=w12+w22+⋯+wn2\| \boldsymbol w \| = \sqrt{w^2_1+w^2_2+\cdots +w^2_n} ∥w∥=w12​+w22​+⋯+wn2​​在上图中记为长度aaa。超平面上的所有向量在法向量方向上的投影长度都相同，记为lll，即原点到此超平面的距离，所以由内积定义可知，超平面上的任一向量x=[x1x2⋯xn]T\boldsymbol x = [\begin{matrix} x_1 & x_2 & \cdots & x_n\end{matrix}]^Tx=[x1​​x2​​⋯​xn​​]T与此法向量w\boldsymbol ww的内积绝对值就是其在法向量方向投影与法向量模长的乘积，即

∣wTx∣=∣w1x1+w2x2+⋯+wnxn∣=al|\boldsymbol w^T \boldsymbol x | = |w_1x_1+w_2x_2+ \cdots + w_nx_n|=al ∣wTx∣=∣w1​x1​+w2​x2​+⋯+wn​xn​∣=al由于一旦法向量和超平面固定，aaa和lll就都是固定值，所以可将上式中alalal改写成−b-b−b，并赋予b>0b>0b>0的含义是向量x\boldsymbol xx投影方向与超平面在法向量反方向，于是就得到了超平面方程：

wTx+b=0\boldsymbol w^T \boldsymbol x + b = 0 wTx+b=0并称其中bbb为超平面的截距。

对于超平面外一向量y=[y1y2⋯yn]T\boldsymbol y = [\begin{matrix} y_1 & y_2 & \cdots & y_n\end{matrix}]^Ty=[y1​​y2​​⋯​yn​​]T，同理可知该向量与超平面法向量的内积绝对值为：

∣wTy∣=∣w1y1+w2y2+⋯+wnyn∣=ad|\boldsymbol w^T \boldsymbol y| = |w_1y_1+w_2y_2+ \cdots + w_ny_n|=ad ∣wTy∣=∣w1​y1​+w2​y2​+⋯+wn​yn​∣=ad而在上图中可直观看出向量y\boldsymbol yy到此超平面的距离是hhh，且有如下关系：

h=d−l=ad−ala=∣wTy+b∣∥w∥h=d-l={ad-al \over a}={|\boldsymbol w^T \boldsymbol y+b| \over \| \boldsymbol w \|} h=d−l=aad−al​=∥w∥∣wTy+b∣​由此可见，超平面外一向量到超平面的几何距离就等于将该向量代入超平面方程取绝对值后，再除以超平面法向量的模长。

（二）严格线性可分问题

设nnn维向量空间中有mmm个向量x1∼xm\boldsymbol x_1 \sim \boldsymbol x_mx1​∼xm​，每个向量xi\boldsymbol x_ixi​都对应一个标签yiy_iyi​，其中一部分向量的标签值是1，另一部分向量的标签值是-1，从而区分出两组向量。我们所要做的是在nnn维向量空间中构造一个超平面 ，使得在超平面两侧找到的支持向量（即距离超平面几何距离最短的向量）xs\boldsymbol x_sxs​到此超平面的距离尽可能远，由前述超平面方程可知这个距离为：

ds=∣wTxs+b∣∥w∥d_s = {|\boldsymbol w^T \boldsymbol x_s + b| \over \| \boldsymbol w \|} ds​=∥w∥∣wTxs​+b∣​同时要满足同标记的向量只能在此超平面的同一边，即

{wTxi+b>0yi=1wTxi+b<0yi=−1\begin{cases} \boldsymbol w^T \boldsymbol x_i + b > 0 & y_i=1 \\ \boldsymbol w^T \boldsymbol x_i + b < 0 & y_i=-1 \end{cases} {wTxi​+b>0wTxi​+b<0​yi​=1yi​=−1​满足上述条件的超平面称为决策面。

二、目标解析

对于上述约束条件，当yi=1y_i=1yi​=1时有wTxi+b>0\boldsymbol w^T \boldsymbol x_i + b > 0wTxi​+b>0，等价于必然存在α>0\alpha>0α>0使得wTxi+b>α\boldsymbol w^T \boldsymbol x_i + b > \alphawTxi​+b>α；同理，必然存在−α-\alpha−α使得当yi=−1y_i=-1yi​=−1时，wTxi+b<−α\boldsymbol w^T \boldsymbol x_i + b < -\alphawTxi​+b<−α；而由超平面方程可知，w\boldsymbol ww和bbb同比例变化后，方程所描述的仍是同一个超平面，所以上述条件等价于不等式两边同除α\alphaα的结果，即

{wTxi+b>1yi=1wTxi+b<−1yi=−1\begin{cases} \boldsymbol w^T \boldsymbol x_i + b > 1 & y_i=1 \\ \boldsymbol w^T \boldsymbol x_i + b < -1 & y_i=-1 \end{cases} {wTxi​+b>1wTxi​+b<−1​yi​=1yi​=−1​由于所追求的决策面到其两侧支持向量（xs+\boldsymbol x_{s+}xs+​和xs−\boldsymbol x_{s-}xs−​）的距离相同，可以令上述α\alphaα就等于支持向量代入决策面表达式后的结果，即

{wTxs++b=αys+=1wTxs−+b=−αys−=−1\begin{cases} \boldsymbol w^T \boldsymbol x_{s+} + b = \alpha & y_{s+}=1 \\ \boldsymbol w^T \boldsymbol x_{s-} + b = -\alpha & y_{s-}=-1 \end{cases} {wTxs+​+b=αwTxs−​+b=−α​ys+​=1ys−​=−1​上述条件同样等价于等式两边同除α\alphaα的结果，即

{wTxs++b=1ys+=1wTxs−+b=−1ys−=−1\begin{cases} \boldsymbol w^T \boldsymbol x_{s+} + b = 1 & y_{s+}=1 \\ \boldsymbol w^T \boldsymbol x_{s-} + b = -1 & y_{s-}=-1 \end{cases} {wTxs+​+b=1wTxs−​+b=−1​ys+​=1ys−​=−1​上述各式等号与不等号两侧再同乘标签值，最终约束条件可简化为:

{yi(wTxi+b)>1i≠s+或s−yi(wTxi+b)=1i=s+或s−\begin{cases} y_i(\boldsymbol w^T \boldsymbol x_i + b) > 1 & i \not= s+ 或 s-\\ y_i(\boldsymbol w^T \boldsymbol x_i + b) = 1 & i = s+ 或 s- \end{cases} {yi​(wTxi​+b)>1yi​(wTxi​+b)=1​i​=s+或s−i=s+或s−​因此，要最大化的目标函数可化为：

ds=∣wTxs+b∣∥w∥=1∥w∥d_s = {|\boldsymbol w^T \boldsymbol x_s + b| \over \| \boldsymbol w \|} = {1 \over \| \boldsymbol w \|} ds​=∥w∥∣wTxs​+b∣​=∥w∥1​综上所述，原始目标可以具体描述为如下求条件最值问题：

{max⁡w1∥w∥s.t.yi(wTxi+b)≥1i=1,2,⋯,m\begin{cases} \max \limits_{\boldsymbol w} {1 \over \| \boldsymbol w \|} \\ s.t. \quad y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1\\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎨⎪⎧​wmax​∥w∥1​s.t.yi​(wTxi​+b)≥1i=1,2,⋯,m​其中s.t.s.t.s.t.是satisfy to的缩写，即同时满足于的意思。

三、兼容软间隔

在实际问题中，总有可能偶尔出现异常数据，比如个别正标向量跑到了负标向量堆里去了，如果严格分类会造成一部分向量与决策面很近，甚至根本无法分界。如下图所示：

为了容忍个别异常数据，提升模型鲁棒性，实现下图分类效果：

需要对原条件最值问题做调整如下：

{max⁡w{1∥w∥−C∑i=1mξi}s.t.yi(wTxi+b)≥1−ξiξi≥0i=1,2,⋯,m\begin{cases} \max \limits_{\boldsymbol w} \{{1 \over \| \boldsymbol w \|}-C \sum \limits_{i=1}^m \xi_i\} \\ s.t. \quad y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1-\xi_i\\ \quad \quad \,\,\, \xi_i \geq 0 \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​wmax​{∥w∥1​−Ci=1∑m​ξi​}s.t.yi​(wTxi​+b)≥1−ξi​ξi​≥0i=1,2,⋯,m​其中新增了两个大于0的变量：

松弛变量ξi\xi_iξi​。当ξi=0\xi_i=0ξi​=0时，说明第iii个向量要严格遵守分类要求；当0<ξi<10 < \xi_i<10<ξi​<1时，说明允许第iii个向量离决策面的距离比支持向量还近；当ξi≥1\xi_i \geq 1ξi​≥1时，表示允许第iii个向量可以处在决策面上甚至在决策面另一边。惩罚系数CCC。如果松弛变量的取值不被约束，则相当于此优化问题不再有约束条件了，这样显然不行，因此需要使每个松弛变量都付出代价，即对优化问题产生负面影响，这个影响度就是惩罚系数，该系数越大则越不能容忍例外向量，越小则约能容忍。故而该系数是超参数，即需要人为设定，而非模型训练出来。

四、兼容非线性分割

当出现类似下图的分类问题时，无论如何容忍异常向量，也无法做到在平面上用一条直线分割两类向量。

对于上图中的向量分布，实际上需要一条环线才能有效分割两类向量。以平面上二次曲线为例，其方程形式为：

w0+w1x1+w2x2+w3x12+w4x22+w5x1x2=0w_0+w_1x_1+w_2x_2+w_3x^2_1+w_4x^2_2+w_5x_1x_2=0 w0​+w1​x1​+w2​x2​+w3​x12​+w4​x22​+w5​x1​x2​=0就相当于将原二维向量映射到五维空间后的线性方程，即

wTϕ(x)+b=0\boldsymbol w^T \phi(\boldsymbol x)+b=0 wTϕ(x)+b=0其中

ϕ(x)=[x1x2x12x22x1x2],w=[w1w2w3w4w5],b=w0\phi (\boldsymbol x) = \left [ \begin{matrix} x_1 \\ x_2 \\ x^2_1 \\ x^2_2 \\ x_1x_2 \end{matrix} \right], \quad \boldsymbol w = \left [ \begin{matrix} w_1 \\ w_2 \\ w_3 \\ w_4 \\ w_5 \end{matrix} \right], \quad b=w_0 ϕ(x)=⎣⎢⎢⎢⎢⎡​x1​x2​x12​x22​x1​x2​​⎦⎥⎥⎥⎥⎤​,w=⎣⎢⎢⎢⎢⎡​w1​w2​w3​w4​w5​​⎦⎥⎥⎥⎥⎤​,b=w0​由此可见，在低维空间无法线性分割的两组向量可以映射到高维空间实现线性分割。为兼容非线性分割问题，需要再次对上述条件最值问题做调整，得到最终目标如下：

{max⁡w{1∥w∥−C∑i=1mξi}s.t.yi[wTϕ(xi)+b]≥1−ξiξi≥0i=1,2,⋯,m\begin{cases} \max \limits_{\boldsymbol w} \{{1 \over \| \boldsymbol w \|}-C \sum \limits_{i=1}^m \xi_i\} \\ s.t. \quad y_i[\boldsymbol w^T \phi (\boldsymbol x_i) + b] \geq 1-\xi_i\\ \quad \quad \,\,\, \xi_i \geq 0 \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​wmax​{∥w∥1​−Ci=1∑m​ξi​}s.t.yi​[wTϕ(xi​)+b]≥1−ξi​ξi​≥0i=1,2,⋯,m​至此，终于将既兼容软间隔，又兼容非线性分割的SVM模型描述成了条件最值问题。后面就是如何求解此条件最值问题了。

五、拉格朗日乘数法

为便于构造拉格朗日函数后计算梯度，可将原条件最值问题等价为：

{min⁡w{∥w∥22+C∑i=1mξi}s.t.1−ξi−yi[wTϕ(xi)+b]≤0ξi≥0i=1,2,⋯,m\begin{cases} \min \limits_{\boldsymbol w} \{{ \| \boldsymbol w \|^2 \over 2}+C \sum \limits_{i=1}^m \xi_i\} \\ s.t. \quad 1-\xi_i - y_i[\boldsymbol w^T \phi (\boldsymbol x_i) + b] \leq 0 \\ \quad \quad \,\,\, \xi_i \geq 0 \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​wmin​{2∥w∥2​+Ci=1∑m​ξi​}s.t.1−ξi​−yi​[wTϕ(xi​)+b]≤0ξi​≥0i=1,2,⋯,m​拉格朗日函数的作用是可以把所有约束条件囊括进一个多元函数，通过求该多元函数的最值来得到想要的条件最值。具体原理证明不在此文中详述，我们只需知道所谓拉格朗日函数就是将所有以小于零作为条件的各条件函数分别乘上一个系数，再加上要求最小值的目标函数后组成的多项式。针对上述条件最值问题，可构造拉格朗日函数如下：

L(w,ξ,b,a,μ)=12∥w∥2+C∑i=1mξi+∑i=1mai{1−ξi−yi[wTϕ(xi)+b]}−∑i=1mμiξiL(\boldsymbol w, \boldsymbol \xi, b, \boldsymbol a, \boldsymbol \mu) = {1 \over 2}\| \boldsymbol w \|^2 + C \sum \limits_{i=1}^m \xi_i + \sum \limits_{i=1}^m a_i\{ 1-\xi_i-y_i[\boldsymbol w^T\phi(\boldsymbol x_i)+b] \} - \sum \limits_{i=1}^m \mu_i \xi_i L(w,ξ,b,a,μ)=21​∥w∥2+Ci=1∑m​ξi​+i=1∑m​ai​{1−ξi​−yi​[wTϕ(xi​)+b]}−i=1∑m​μi​ξi​由于优化的目标是求条件最小值，而约束条件又是一系列不等式，如果直接对上述LLL函数取梯度为零点，则相当于在约束边界上找最小值，但我们想要的效果是最小值只要在约束区域内就行，不一定只能在边界上。也就是说通过LLL函数中(w,ξ,b)(\boldsymbol w, \boldsymbol \xi, b)(w,ξ,b)变量（即除去约束条件前的系数变量）求得的最小值位置可能落在约束区域内也可能落在外面，对于这两种情况，我们想要的效果是：

如果恰好就落在哪几个约束区域内，则这几个约束条件就可以无视了，即其前面的条件系数就可以为0；如果落在哪几个约束区域外，则这几个约束条件就要起作用，把取值点限定在只能在该条件边界上。

要实现上述目的，我们只需要先在限定条件系数（即a\boldsymbol aa和μ\boldsymbol \muμ）大于0的同时由这些条件系数求LLL函数最大值后，再通过非条件系数变量(w,ξ,b)(\boldsymbol w, \boldsymbol \xi, b)(w,ξ,b)求LLL函数最小值即可。因为约束条件都是约束小于0的，所以

一旦某个约束条件被满足了，则该条件项就变成负数了，此时要想让LLL函数值尽可能大，则该条件项前的条件系数就只能变成0，这样第一个想要的效果就满足了；一旦某个约束条件没被满足，则该条件项就变成正数了，此时要想让LLL函数值尽可能大，则该条件项前的条件系数就会尽可能大，甚至到正无穷。而后续还会对LLL函数取最小值，这样就使得这些项前系数取到正无穷的条件项只能取0，即把取值点限定在该条件边界上，从而也满足了第二个想要的效果。

因此，最终的目标可描述为：

min⁡w,ξ,bmax⁡a≥0,μ≥0L(w,ξ,b,a,μ)\min \limits_{\boldsymbol w, \boldsymbol \xi, b} \max \limits_{\boldsymbol a \geq 0, \boldsymbol \mu \geq 0} L(\boldsymbol w, \boldsymbol \xi, b, \boldsymbol a, \boldsymbol \mu) w,ξ,bmin​a≥0,μ≥0max​L(w,ξ,b,a,μ)由于取最大结果中的最小值和取最小结果中的最大值效果相同（即拉格朗日对偶性），因此上述问题又等价于:

max⁡a≥0,μ≥0min⁡w,ξ,bL(w,ξ,b,a,μ)\max \limits_{\boldsymbol a \geq 0, \boldsymbol \mu \geq 0} \min \limits_{\boldsymbol w, \boldsymbol \xi, b} L(\boldsymbol w, \boldsymbol \xi, b, \boldsymbol a, \boldsymbol \mu) a≥0,μ≥0max​w,ξ,bmin​L(w,ξ,b,a,μ)第一步先取最小值，通过使LLL函数分别对各自变量的偏导为零可得，即

∇w,ξ,bL={∂L∂w=w−∑i=1maiyiϕ(xi)=0∂L∂b=−∑i=1maiyi=0∂L∂ξi=C−ai−μi=0\nabla_{\boldsymbol w, \boldsymbol \xi, b}L= \begin{cases} {\partial L \over \partial \boldsymbol w} = \boldsymbol w - \sum \limits_{i=1}^m a_iy_i \phi (\boldsymbol x_i) = 0 \\ {\partial L \over \partial b} = - \sum \limits_{i=1}^m a_iy_i = 0 \\ {\partial L \over \partial \xi_i} = C - a_i - \mu_i = 0 \end{cases} ∇w,ξ,b​L=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​∂w∂L​=w−i=1∑m​ai​yi​ϕ(xi​)=0∂b∂L​=−i=1∑m​ai​yi​=0∂ξi​∂L​=C−ai​−μi​=0​将上述结论代入第一步优化目标得

min⁡w,ξ,bL(w,ξ,b,a,μ)=12∥w∥2+C∑i=1mξi+∑i=1mai{1−ξi−yi[wTϕ(xi)+b]}−∑i=1mμiξi=12wTw+C∑i=1mξi+∑i=1mai−∑i=1maiξi−wTw−b∑i=1maiyi−C∑i=1mξi+∑i=1maiξi=∑i=1mai−12wTw=∑i=1mai−12∑i=1m∑j=1maiajyiyjϕ(xi)Tϕ(xj)\begin{aligned} \min \limits_{\boldsymbol w, \boldsymbol \xi, b} L(\boldsymbol w, \boldsymbol \xi, b, \boldsymbol a, \boldsymbol \mu) &= {1 \over 2}\| \boldsymbol w \|^2 + C \sum \limits_{i=1}^m \xi_i + \sum \limits_{i=1}^m a_i\{ 1-\xi_i-y_i[\boldsymbol w^T\phi(\boldsymbol x_i)+b] \} - \sum \limits_{i=1}^m \mu_i \xi_i \\ &= {1 \over 2} \boldsymbol w^T \boldsymbol w + C \sum \limits_{i=1}^m \xi_i + \sum \limits_{i=1}^m a_i - \sum \limits_{i=1}^m a_i \xi_i - \boldsymbol w^T \boldsymbol w - b \sum \limits_{i=1}^m a_i y_i - C\sum \limits_{i=1}^m \xi_i +\sum \limits_{i=1}^m a_i \xi_i \\ &= \sum \limits_{i=1}^m a_i - {1 \over 2}\boldsymbol w^T \boldsymbol w \\ &= \sum \limits_{i=1}^m a_i - {1 \over 2}\sum \limits_{i=1}^m \sum \limits_{j=1}^m a_i a_j y_i y_j \phi(\boldsymbol x_i)^T \phi(\boldsymbol x_j) \end{aligned}w,ξ,bmin​L(w,ξ,b,a,μ)​=21​∥w∥2+Ci=1∑m​ξi​+i=1∑m​ai​{1−ξi​−yi​[wTϕ(xi​)+b]}−i=1∑m​μi​ξi​=21​wTw+Ci=1∑m​ξi​+i=1∑m​ai​−i=1∑m​ai​ξi​−wTw−bi=1∑m​ai​yi​−Ci=1∑m​ξi​+i=1∑m​ai​ξi​=i=1∑m​ai​−21​wTw=i=1∑m​ai​−21​i=1∑m​j=1∑m​ai​aj​yi​yj​ϕ(xi​)Tϕ(xj​)​第二步就是进一步求解max⁡a≥0,μ≥0min⁡w,ξ,bL(w,ξ,b,a,μ)\max \limits_{\boldsymbol a \geq 0, \boldsymbol \mu \geq 0} \min \limits_{\boldsymbol w, \boldsymbol \xi, b} L(\boldsymbol w, \boldsymbol \xi, b, \boldsymbol a, \boldsymbol \mu)a≥0,μ≥0max​w,ξ,bmin​L(w,ξ,b,a,μ)，结合上述条件，原条件最值问题可化为如下优化问题：

{min⁡a{12∑i=1m∑j=1maiajyiyjϕ(xi)Tϕ(xj)−∑i=1mai}s.t.∑i=1maiyi=0ai,μi≥0C−ai−μi=0i=1,2,⋯,m\begin{cases} \min \limits_{\boldsymbol a} \{{1 \over 2}\sum \limits_{i=1}^m \sum \limits_{j=1}^m a_i a_j y_i y_j \phi(\boldsymbol x_i)^T \phi(\boldsymbol x_j) - \sum \limits_{i=1}^m a_i\} \\ s.t. \quad \sum \limits_{i=1}^m a_iy_i = 0 \\ \quad \quad \,\,\, a_i,\mu_i \geq 0 \\ \quad \quad \,\,\, C - a_i - \mu_i = 0 \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​amin​{21​i=1∑m​j=1∑m​ai​aj​yi​yj​ϕ(xi​)Tϕ(xj​)−i=1∑m​ai​}s.t.i=1∑m​ai​yi​=0ai​,μi​≥0C−ai​−μi​=0i=1,2,⋯,m​由于

{μi≥0μi=C−ai⇒ai≤C\begin{cases} \mu_i \geq 0 \\ \mu_i = C - a_i \end{cases} \Rightarrow a_i \leq C {μi​≥0μi​=C−ai​​⇒ai​≤C从而得到最终条件最值问题：

{min⁡a{12∑i=1m∑j=1maiajyiyjϕ(xi)Tϕ(xj)−∑i=1mai}s.t.∑i=1maiyi=00≤ai≤Ci=1,2,⋯,m\begin{cases} \min \limits_{\boldsymbol a} \{{1 \over 2}\sum \limits_{i=1}^m \sum \limits_{j=1}^m a_i a_j y_i y_j \phi(\boldsymbol x_i)^T \phi(\boldsymbol x_j) - \sum \limits_{i=1}^m a_i\} \\ s.t. \quad \sum \limits_{i=1}^m a_iy_i = 0 \\ \quad \quad \,\,\, 0 \leq a_i \leq C \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​amin​{21​i=1∑m​j=1∑m​ai​aj​yi​yj​ϕ(xi​)Tϕ(xj​)−i=1∑m​ai​}s.t.i=1∑m​ai​yi​=00≤ai​≤Ci=1,2,⋯,m​上述问题涉及映射到高维空间的向量内积，即ϕ(xi)Tϕ(xj)\phi(\boldsymbol x_i)^T \phi(\boldsymbol x_j)ϕ(xi​)Tϕ(xj​)的计算问题，需要引入核函数。

六、核函数

在第四章的例子中，对二维向量做二次曲线分割就需要映射到五维空间，如果三维甚至更高维度向量做高次超曲面分割，所要升的维度将很高，从而造成极大的计算量。因此需找到一种方法，使得两个低维向量无需向高维空间做映射，只需在原维度空间做简单运算，所得结果就等于其映射到高维空间后的内积。这种方法就是核函数，将同是低维度的两个向量代入核函数，所得结果就是其映射到高维空间后的内积，即

ϕ(xi)Tϕ(xj)=κ(xi,xj)\phi(\boldsymbol x_i)^T \phi(\boldsymbol x_j) = \kappa(\boldsymbol x_i, \boldsymbol x_j) ϕ(xi​)Tϕ(xj​)=κ(xi​,xj​)常用的核函数有以下几种：

1、线性核。就表示原空间内积，适用于线性可分问题。

κ(xi,xj)=xiTxj\kappa(\boldsymbol x_i, \boldsymbol x_j) = \boldsymbol x_i^T \boldsymbol x_j κ(xi​,xj​)=xiT​xj​2、高斯核。适用于没有先验经验的非线性分类。其中σ\sigmaσ越小，使得映射的维度越高。

κ(xi,xj)=exp⁡(−∥xi−xj∥22σ2)\kappa(\boldsymbol x_i, \boldsymbol x_j) = \exp \left( -{\|\boldsymbol x_i - \boldsymbol x_j \|^2 \over 2\sigma^2} \right ) κ(xi​,xj​)=exp(−2σ2∥xi​−xj​∥2​)3、多项式核。适用于没有先验经验的分类。其中ddd越越小，使得映射的维度越高。

κ(xi,xj)=(axiTxj+c)d\kappa(\boldsymbol x_i, \boldsymbol x_j) =( a \boldsymbol x^T_i \boldsymbol x_j +c)^d κ(xi​,xj​)=(axiT​xj​+c)d4、Sigmoid核。此时SVM实现的就是一种多层感知器神经网络。

κ(xi,xj)=tanh⁡(axiTxj+c)\kappa(\boldsymbol x_i, \boldsymbol x_j) =\tanh( a \boldsymbol x^T_i \boldsymbol x_j +c) κ(xi​,xj​)=tanh(axiT​xj​+c)

大多数非线性分割使用高斯核都能很好处理。

七、SMO算法

通过核函数解决了高维向量内积计算问题后，就剩下如何求解出满足约束条件的一组最优条件系数a\boldsymbol aa这最后一个问题了。由于条件系数aia_iai​的个数与待分割向量一样多，难以直接计算梯度最低点，故而采用SMO算法，即序列最小优化（Sequential Minimal Optimization）算法。其思路是每次选择两个待优化条件系数而固定其他系数，由约束条件又可以化去一个系数，从而化成一个一元二次函数的求最值问题。对于已满足约束条件的最优解，再使用SMO算法做迭代优化时，是不会改变其结果的，所以当通过两两优化的方式把所有条件系数遍历几遍，各系数已不再有可优化空间后，即得到了最优结果。

（一）无约束最小化

首先选定a1a_1a1​和a2a_2a2​做优化，其他系数由于都先被视为常量，所以可以都先无视；并用核函数代替映射到高维空间的向量内积，则前述条件最值问题可化为：

{min⁡a{12a12κ(x1,x1)+12a22κ(x2,x2)+y1y2a1a2κ(x1,x2)+y1a1∑i=3myiaiκ(xi,x1)+y2a2∑i=3myiaiκ(xi,x2)−(a1+a2)}s.t.a1y1+a2y2=−∑i=3maiyi=ς0≤ai≤Ci=1,2,⋯,m\begin{cases} \min \limits_{\boldsymbol a} \{{1 \over 2}a_1^2 \kappa(\boldsymbol x_1, \boldsymbol x_1) + {1 \over 2}a_2^2 \kappa(\boldsymbol x_2, \boldsymbol x_2) + y_1y_2a_1a_2\kappa(\boldsymbol x_1, \boldsymbol x_2) + y_1a_1\sum \limits_{i=3}^m y_i a_i \kappa(\boldsymbol x_i, \boldsymbol x_1) + y_2a_2\sum \limits_{i=3}^m y_i a_i \kappa(\boldsymbol x_i, \boldsymbol x_2) - (a_1+a_2) \} \\ s.t. \quad a_1y_1 + a_2y_2 = - \sum \limits_{i=3}^m a_i y_i = \varsigma\\ \quad \quad \,\,\, 0 \leq a_i \leq C \\ \quad \quad \,\,\, i = 1,2, \cdots,m \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​amin​{21​a12​κ(x1​,x1​)+21​a22​κ(x2​,x2​)+y1​y2​a1​a2​κ(x1​,x2​)+y1​a1​i=3∑m​yi​ai​κ(xi​,x1​)+y2​a2​i=3∑m​yi​ai​κ(xi​,x2​)−(a1​+a2​)}s.t.a1​y1​+a2​y2​=−i=3∑m​ai​yi​=ς0≤ai​≤Ci=1,2,⋯,m​其中ς\varsigmaς为常数；κ(xi,xj)\kappa(\boldsymbol x_i, \boldsymbol x_j)κ(xi​,xj​)可简记为κij\kappa_{ij}κij​。由此，上述要最小化的目标函数可记为：

J=12a12κ11+12a22κ22+y1y2a1a2κ12+y1a1∑i=3myiaiκi1+y2a2∑i=3myiaiκi2−(a1+a2)J={1 \over 2}a_1^2 \kappa_{11} + {1 \over 2}a_2^2 \kappa_{22} + y_1y_2a_1a_2 \kappa_{12} + y_1a_1\sum \limits_{i=3}^m y_i a_i \kappa_{i1} + y_2a_2\sum \limits_{i=3}^m y_i a_i \kappa_{i2} - (a_1+a_2) J=21​a12​κ11​+21​a22​κ22​+y1​y2​a1​a2​κ12​+y1​a1​i=3∑m​yi​ai​κi1​+y2​a2​i=3∑m​yi​ai​κi2​−(a1​+a2​)由于a1y1+a2y2=ςa_1y_1 + a_2y_2 = \varsigmaa1​y1​+a2​y2​=ς，且yiy_iyi​只可能取±1\pm 1±1，可得a1=y1(ς−a2y2)a_1=y_1(\varsigma - a_2y_2)a1​=y1​(ς−a2​y2​)，代入上式对a2a_2a2​求偏导为：

∂J∂a2=(κ11+κ22−2κ12)a2+(y2κ12−y2κ11)ς−y2∑i=3myiaiκi1+y2∑i=3myiaiκi2+y1y2−1{\partial J \over \partial a_2}= (\kappa_{11} + \kappa_{22} - 2 \kappa_{12} ) a_2+ (y_2 \kappa_{12} - y_2 \kappa_{11})\varsigma - y_2\sum \limits_{i=3}^m y_i a_i \kappa_{i1} + y_2\sum \limits_{i=3}^m y_i a_i \kappa_{i2} + y_1y_2 - 1 ∂a2​∂J​=(κ11​+κ22​−2κ12​)a2​+(y2​κ12​−y2​κ11​)ς−y2​i=3∑m​yi​ai​κi1​+y2​i=3∑m​yi​ai​κi2​+y1​y2​−1令上式等于0，同时将待优化更新的a2a_2a2​记为a2′a_2'a2′​，区别于其他仍未更新的旧值aia_iai​，有：

(κ11+κ22−2κ12)a2′=(y2κ11−y2κ12)ς+y2∑i=3myiaiκi1−y2∑i=3myiaiκi2−y1y2+1(\kappa_{11} + \kappa_{22} - 2 \kappa_{12} ) a_2' = (y_2 \kappa_{11} - y_2 \kappa_{12})\varsigma +y_2\sum \limits_{i=3}^m y_i a_i \kappa_{i1} -y_2\sum \limits_{i=3}^m y_i a_i \kappa_{i2} -y_1y_2 + 1 (κ11​+κ22​−2κ12​)a2′​=(y2​κ11​−y2​κ12​)ς+y2​i=3∑m​yi​ai​κi1​−y2​i=3∑m​yi​ai​κi2​−y1​y2​+1由于

∑i=1myiaiκij=a1y1κ1j+a2y2κ2j+∑i=3myiaiκij\sum \limits_{i=1}^m y_i a_i \kappa_{ij} =a_1y_1\kappa_{1j} +a_2y_2\kappa_{2j}+\sum \limits_{i=3}^m y_i a_i \kappa_{ij} i=1∑m​yi​ai​κij​=a1​y1​κ1j​+a2​y2​κ2j​+i=3∑m​yi​ai​κij​代入上式可得：

(κ11+κ22−2κ12)a2′=(y2κ11−y2κ12)ς+y2(∑i=1myiaiκi1−a1y1κ11−a2y2κ21)−y2(∑i=1myiaiκi2−a1y1κ12−a2y2κ22)−y1y2+1=(y2κ11−y2κ12)ς+y2∑i=1myiaiκi1−y1y2a1κ11−a2κ21−y2∑i=1myiaiκi2+y1y2a1κ12+a2κ22−y1y2+1\begin{aligned} (\kappa_{11} + \kappa_{22} - 2 \kappa_{12} ) a_2' &= (y_2 \kappa_{11} - y_2 \kappa_{12})\varsigma +y_2(\sum \limits_{i=1}^m y_i a_i \kappa_{i1} -a_1y_1\kappa_{11} - a_2y_2\kappa_{21}) -y_2(\sum \limits_{i=1}^m y_i a_i \kappa_{i2} -a_1y_1\kappa_{12} - a_2y_2\kappa_{22}) -y_1y_2 + 1 \\ &=(y_2 \kappa_{11} - y_2 \kappa_{12})\varsigma +y_2\sum \limits_{i=1}^m y_i a_i \kappa_{i1} - y_1y_2a_1\kappa_{11} - a_2\kappa_{21} -y_2\sum \limits_{i=1}^m y_i a_i \kappa_{i2} +y_1y_2a_1\kappa_{12} + a_2\kappa_{22} -y_1y_2 + 1 \end{aligned}(κ11​+κ22​−2κ12​)a2′​​=(y2​κ11​−y2​κ12​)ς+y2​(i=1∑m​yi​ai​κi1​−a1​y1​κ11​−a2​y2​κ21​)−y2​(i=1∑m​yi​ai​κi2​−a1​y1​κ12​−a2​y2​κ22​)−y1​y2​+1=(y2​κ11​−y2​κ12​)ς+y2​i=1∑m​yi​ai​κi1​−y1​y2​a1​κ11​−a2​κ21​−y2​i=1∑m​yi​ai​κi2​+y1​y2​a1​κ12​+a2​κ22​−y1​y2​+1​为简化表达式，记

Ej=∑i=1myiaiκij−yjE_j = \sum \limits_{i=1}^m y_i a_i \kappa_{ij} - y_j Ej​=i=1∑m​yi​ai​κij​−yj​即

E=[E1E2⋮Em]=[∑i=1myiaiκi1−y1∑i=1myiaiκi2−y2⋮∑i=1myiaiκim−ym]=[κ11κ12⋯κ1mκ21κ22⋯κ2m⋮⋮⋱⋮κm1κm2⋯κmm][y1a1y2a2⋮ymam]−[y1y2⋮ym]\boldsymbol E = \left [ \begin{matrix} E_1\\ E_2 \\ \vdots \\ E_m \end{matrix} \right] = \left [ \begin{matrix} \sum \limits_{i=1}^m y_i a_i \kappa_{i1} - y_1 \\ \sum \limits_{i=1}^m y_i a_i \kappa_{i2} - y_2 \\ \vdots \\ \sum \limits_{i=1}^m y_i a_i \kappa_{im} - y_m \end{matrix} \right] = \left [ \begin{matrix} \kappa_{11} & \kappa_{12} & \cdots & \kappa_{1m} \\ \kappa_{21} & \kappa_{22} & \cdots & \kappa_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ \kappa_{m1} & \kappa_{m2} & \cdots & \kappa_{mm} \\ \end{matrix} \right] \left [ \begin{matrix} y_1 a_1 \\ y_2 a_2 \\ \vdots \\ y_m a_m \end{matrix} \right] - \left [ \begin{matrix} y_1\\ y_2 \\ \vdots \\ y_m \end{matrix} \right] E=⎣⎢⎢⎢⎡​E1​E2​⋮Em​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​i=1∑m​yi​ai​κi1​−y1​i=1∑m​yi​ai​κi2​−y2​⋮i=1∑m​yi​ai​κim​−ym​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​κ11​κ21​⋮κm1​​κ12​κ22​⋮κm2​​⋯⋯⋱⋯​κ1m​κ2m​⋮κmm​​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​y1​a1​y2​a2​⋮ym​am​​⎦⎥⎥⎥⎤​−⎣⎢⎢⎢⎡​y1​y2​⋮ym​​⎦⎥⎥⎥⎤​则有：

y2(E1−E2）=y2∑i=1myiaiκi1−y1y2−y2∑i=1myiaiκi2+1y_2(E_1-E2） = y_2\sum \limits_{i=1}^m y_i a_i \kappa_{i1} - y_1y_2 - y_2\sum \limits_{i=1}^m y_i a_i \kappa_{i2} + 1 y2​(E1​−E2）=y2​i=1∑m​yi​ai​κi1​−y1​y2​−y2​i=1∑m​yi​ai​κi2​+1代入前式可化简为：

(κ11+κ22−2κ12)a2′=(y2κ11−y2κ12)ς−y1y2a1κ11−a2κ21+y1y2a1κ12+a2κ22+y2(E1−E2)(\kappa_{11} + \kappa_{22} - 2 \kappa_{12} ) a_2' = (y_2 \kappa_{11} - y_2 \kappa_{12})\varsigma -y_1y_2a_1\kappa_{11} - a_2\kappa_{21} +y_1y_2a_1\kappa_{12} + a_2\kappa_{22} +y_2(E_1-E_2) (κ11​+κ22​−2κ12​)a2′​=(y2​κ11​−y2​κ12​)ς−y1​y2​a1​κ11​−a2​κ21​+y1​y2​a1​κ12​+a2​κ22​+y2​(E1​−E2​)再将a1=y1(ς−a2y2)a_1=y_1(\varsigma - a_2y_2)a1​=y1​(ς−a2​y2​)代入上式可得：

(κ11+κ22−2κ12)a2′=y2κ11ς−y2κ12ς−y2(ς−a2y2)κ11−a2κ21+y2(ς−a2y2)κ12+a2κ22+y2(E1−E2)=a2κ11−a2κ21−a2κ12+a2κ22+y2(E1−E2)=(κ11+κ22−2κ12)a2+y2(E1−E2)\begin{aligned} (\kappa_{11} + \kappa_{22} - 2 \kappa_{12} ) a_2' &= y_2 \kappa_{11} \varsigma - y_2 \kappa_{12}\varsigma -y_2(\varsigma - a_2y_2)\kappa_{11} - a_2\kappa_{21} +y_2(\varsigma - a_2y_2)\kappa_{12} + a_2\kappa_{22} +y_2(E_1-E_2) \\ &=a_2 \kappa_{11} - a_2 \kappa_{21} - a_2 \kappa_{12} + a_2 \kappa_{22} + y_2(E_1-E_2) \\ &=(\kappa_{11} + \kappa_{22} - 2 \kappa_{12})a_2 + y_2(E_1-E_2) \end{aligned}(κ11​+κ22​−2κ12​)a2′​​=y2​κ11​ς−y2​κ12​ς−y2​(ς−a2​y2​)κ11​−a2​κ21​+y2​(ς−a2​y2​)κ12​+a2​κ22​+y2​(E1​−E2​)=a2​κ11​−a2​κ21​−a2​κ12​+a2​κ22​+y2​(E1​−E2​)=(κ11​+κ22​−2κ12​)a2​+y2​(E1​−E2​)​因此：

a2′=a2+y2E1−E2κ11+κ22−2κ12a_2' = a_2 + y_2{E_1-E_2 \over \kappa_{11} + \kappa_{22} - 2 \kappa_{12}} a2′​=a2​+y2​κ11​+κ22​−2κ12​E1​−E2​​上式右侧a2a_2a2​是未更新的旧值。此时得到的a2′a_2'a2′​还并没有考虑约束条件0≤ai≤C0 \leq a_i \leq C0≤ai​≤C ，因此还有待修剪。为方便区分，最终修剪后的各条件系数记为a^i\hat a_ia^i​

（二）根据约束条件做修剪

结合约束条件0≤ai≤C0 \leq a_i \leq C0≤ai​≤C及a1y1+a2y2=ςa_1y_1 + a_2y_2 = \varsigmaa1​y1​+a2​y2​=ς，可确定(a^1,a^2)(\hat a_1, \hat a_2)(a^1​,a^2​)的坐标只能落在下图方框中：

由y1y_1y1​和y2y_2y2​的不同取值可分四种情况讨论：

当y1=y2=1y_1=y_2=1y1​=y2​=1时，有a^2=−a^1+ς\hat a_2 = -\hat a_1+\varsigmaa^2​=−a^1​+ς。此时ς≥0\varsigma \geq 0ς≥0，如上面左图所示，区分ς\varsigmaς与CCC的位置关系做分类讨论如下：

（1）当0≤ς≤C0 \leq \varsigma \leq C0≤ς≤C时，a^1∈[0,ς]\hat a_1 \in [0, \varsigma]a^1​∈[0,ς]，所以a^2∈[0,ς=a1+a2]\hat a_2 \in [0, \varsigma=a_1+a_2]a^2​∈[0,ς=a1​+a2​]。

（2）当C≤ς≤2CC \leq \varsigma \leq 2CC≤ς≤2C时，a^1∈[ς−C,C]\hat a_1 \in [\varsigma-C, C]a^1​∈[ς−C,C]，所以a^2∈[ς−C=a1+a2−C,C]\hat a_2 \in [\varsigma-C=a_1+a_2-C, C]a^2​∈[ς−C=a1​+a2​−C,C]。当y1=y2=−1y_1=y_2=-1y1​=y2​=−1时，有a^2=−a^1−ς\hat a_2 = -\hat a_1-\varsigmaa^2​=−a^1​−ς。此时ς≤0\varsigma \leq 0ς≤0，将此处取负值的ς\varsigmaς换成上述取正值的ς\varsigmaς，可得与上述相同的结果。当y1=1y_1=1y1​=1，y2=−1y_2=-1y2​=−1时，有a^2=a^1−ς\hat a_2 = \hat a_1-\varsigmaa^2​=a^1​−ς。此时−C≤ς≤C-C \leq \varsigma \leq C−C≤ς≤C，如上面右图所示，区分ς\varsigmaς与000的位置关系做分类讨论如下：

（1）当−C≤ς≤0-C \leq \varsigma \leq 0−C≤ς≤0时，a^1∈[0,C+ς]\hat a_1 \in [0, C+\varsigma]a^1​∈[0,C+ς]，所以a^2∈[−ς=a2−a1,C]\hat a_2 \in [-\varsigma=a_2-a_1, C]a^2​∈[−ς=a2​−a1​,C]。

（2）当0≤ς≤C0 \leq \varsigma \leq C0≤ς≤C时，a^1∈[ς,C]\hat a_1 \in [\varsigma, C]a^1​∈[ς,C]，所以a^2∈[0,C−ς=C+a2−a1]\hat a_2 \in [0, C-\varsigma=C+a_2-a_1]a^2​∈[0,C−ς=C+a2​−a1​]。当y1=−1y_1=-1y1​=−1，y2=1y_2=1y2​=1时，有a^2=a^1+ς\hat a_2 = \hat a_1 + \varsigmaa^2​=a^1​+ς。同理可将此处ς\varsigmaς换成−ς-\varsigma−ς，可得与上述相同的结果。

综上可得a^2\hat a_2a^2​的取值上下界为：

{if(y1=y2)L=max⁡{0,a1+a2−C}H=min⁡{a1+a2,C}if(y1≠y2)L=max⁡{a2−a1,0}H=min⁡{C,C+a2−a1}\begin{cases} if(y_1 = y_2) \quad L=\max\{ 0, a_1+a_2-C\} & H=\min\{ a_1+a_2,C\} \\ if(y_1 \not = y_2) \quad L=\max\{ a_2-a_1,0\} & H=\min\{C, C+a_2-a_1\} \end{cases} {if(y1​=y2​)L=max{0,a1​+a2​−C}if(y1​​=y2​)L=max{a2​−a1​,0}​H=min{a1​+a2​,C}H=min{C,C+a2​−a1​}​注意，上述a1a_1a1​和a2a_2a2​都是指本轮更新前的原值。至此可得：

a^2={Ha2′>Ha2′L≤a2′≤HLa2′<Ha^1=y1(ς−a^2y2)=y1[(y1a1+y2a2)−a^2y2]=a1+y1y2(a2−a^2)\begin{aligned} \hat a_2 &= \begin{cases} H & a_2' > H \\ a_2' & L \leq a_2' \leq H \\ L & a_2' < H \\ \end{cases} \\ \hat a_1 &= y_1(\varsigma-\hat a_2 y_2)=y_1[(y_1a_1+y_2a_2)-\hat a_2y_2]=a_1+y_1y_2(a_2-\hat a_2) \end{aligned}a^2​a^1​​=⎩⎪⎨⎪⎧​Ha2′​L​a2′​>HL≤a2′​≤Ha2′​<H​=y1​(ς−a^2​y2​)=y1​[(y1​a1​+y2​a2​)−a^2​y2​]=a1​+y1​y2​(a2​−a^2​)​以此类推，将所有条件系数aia_iai​遍历更新几遍后即得到最终优化后的一组条件系数a^i\hat a_ia^i​。

由第五章可知，在拉格朗日函数中，对于不等式约束条件项，要么该约束条件表达式取0，要么该条件项前的条件系数取0（即所谓KKT条件），即若a^s≠0\hat a_s \not= 0a^s​​=0，则ys(∑i=1maiyiκis+b)≥1−ξsy_s(\sum \limits_{i=1}^m a_iy_i \kappa_{is} + b) \geq 1-\xi_sys​(i=1∑m​ai​yi​κis​+b)≥1−ξs​，也就是说该xs\boldsymbol x_sxs​向量要么是支持向量，即ξs=0\xi_s=0ξs​=0；要么就是被容忍的异常向量，即ξs>0\xi_s>0ξs​>0。对于支持向量，有

b=ys−∑i=1maiyiκis=−Esb = y_s - \sum \limits_{i=1}^m a_iy_i \kappa_{is} = -E_s b=ys​−i=1∑m​ai​yi​κis​=−Es​所以我们只需找到支持向量就可求得决策面截距。同时上式也说明凡支持向量，其对应的偏移量EEE都是相等的，所以我们只需要找到所有条件系数不为零的特殊向量，再从这些特殊向量所对应的偏移量中找到出现次数最多的值，就是支持向量的偏移量，再取负数就得到了决策面的截距b^\hat bb^。最终得到决策面方程为：

∑i=1ma^iyiκ(xi,x)+b^=0\sum \limits_{i=1}^m \hat a_iy_i \kappa(\boldsymbol x_i, \boldsymbol x) + \hat b = 0 i=1∑m​a^i​yi​κ(xi​,x)+b^=0

八、应用示例

（一）SVM函数代码

import numpy as npimport numpy.matlib as matlibimport numpy.random as randomfrom collections import Counter# 定义高斯核函数def gaussKernel(X, Y, sigma=3):'''Parameters----------X : np.matrix一个n行m列矩阵，代表m个n维向量。Y : np.matrix一个n行m列矩阵，代表另外m个n维向量。sigma : float调控参数。越小映射的空间维度约高。默认为3。Returns-------K : np.matrix一个m行1列矩阵，其中k_i = exp[-(||X_i-Y_i||^2)/(2*sigma^2)]。Examples--------X = np.matrix('1, 2, 3; 4, 5, 6').T Y = np.matrix('1, 3, 5; 2, 4, 6').TK = gaussKernel(X, Y, 3)'''K = np.exp(-(np.square(X - Y).T * matlib.ones((np.shape(X)[0], 1))) / (2*sigma**2))return K# 基于SVM对一组做了二元标记的向量给出划分界限的决策平面def svnSimple(dataMatrix, labelVector, C, maxIter, kernel=None, sigma=3):'''Parameters----------dataMatrix : np.matrix一个n行m列矩阵，代表m个n维向量，作为待分类向量。labelVector : np.matrix一个由+1和-1组成的1行m列矩矩阵，作为dataMatrix各列向量的标签。C : float惩罚系数。一个大于等于0的数值，越接近0则对异常向量的容忍度越高。maxIter : int最大遍历次数。对所有条件系数做一次迭代更新为一次遍历。kernel : function核函数。用于在低维空间计算映射到高维空间的向量内积。sigma : float核函数的调控参数。越小映射的空间维度约高。默认为3。Returns-------A : np.matrix一个1行m列矩阵，代表最优条件系数向量。b : float一个1行1列矩阵，代表决策平面的截距。Examples--------X = np.matrix('1, 1; 4, 3; 3, 3').T Y = np.matrix('-1, 1, 1')A, b = svnSimple(X, Y, 100, 10)'''n, m = np.shape(dataMatrix) # 初始化待分向量维度和向量个数A = matlib.zeros((1, m)) # 初始化条件系数向量# 初始化待分向量间内积矩阵if callable(kernel):# 如果给定了核函数则使用核函数做高维内积K = matlib.zeros((m, m)) # 初始化高维内积矩阵for i in range(m):K[:,i] = kernel(dataMatrix[:,i] * matlib.ones((1,m)), dataMatrix, sigma)else:K = dataMatrix.T * dataMatrix# 由SMO算法迭代出所有最优条件系数iterNum = 0# 初始化迭代次数effTraversalNum = 0 # 初始化有效遍历次数while (iterNum < maxIter):alphaPairsChanged = 0 # 初始化条件系数更新次数for i in range(m):# 计算各待分向量到决策面的偏移值（不含b影响）向量E = (K*np.multiply(A, labelVector).T ).T- labelVector# 从其他条件系数中再随机选出一个与当前条件系数作为一对待优化系数j = iwhile (j == i):j = int(random.uniform(0, m))# 计算当前待优化的第二个条件系数a2的待修剪值a2 = A[0,j] + labelVector[0,j]*(E[0,i]-E[0,j])/(K[i,i]+K[j,j]-2*K[i,j])# 计算a2的上下界if labelVector[0,i] == labelVector[0,j]:l = max(0, A[0,i]+A[0,j]-C)h = min(C, A[0,i]+A[0,j])else:l = max(0, A[0,j]-A[0,i])h = min(C, C+A[0,j]-A[0,i])# 修剪条件系数a2if a2 > h:a2 = helif a2 < l:a2 = l# 当a2更新变化太小时便无需更新if (abs(A[0,j] - a2) < 0.00001): continue# 计算条件系数a1a1 = A[0,i] + labelVector[0,i]*labelVector[0,j]*(A[0,j]-a2)# 更新条件系数向量A[0,i] = a1A[0,j] = a2# 统计本此遍历中条件系数更新次数alphaPairsChanged += 1# 统计有效遍历次数if alphaPairsChanged != 0: effTraversalNum += 1# 遍历迭代次数加1iterNum += 1print("共完成有效遍历%d次。" %effTraversalNum)# 通过支持向量计算决策平面截距bspVecIndex = [] # 初始化特殊向量序号集合for k in range(m):if abs(A[0,k]) > 0.01: # 条件系数显著大于零对应的是特殊向量spVecIndex.append(k)spE = E[:, spVecIndex].tolist()[0] # 获取所有特殊向量偏移值roundSpE = [] for n in spE: roundSpE.append(round(n, 4)) # 精确到小数点后4位# 对各偏移量按出现次数由大到小排序listRoundSpE = list(Counter(roundSpE).items())listRoundSpE.sort(key=lambda x:x[1],reverse=True)b = -listRoundSpE[0][0]# 出现次数最多的E对应支持向量return A, b

（二）绘图代码

import numpy as npimport numpy.matlib as matlibimport matplotlib.pyplot as plt# 绘制待分二维向量散点图及其分界线（如有）def showDataSet(dataMatrix, labelVector, A=None, b=None, kernel=None, sigma=None):'''Parameters----------dataMatrix : np.matrix一个2行m列矩阵，代表m个2维向量，作为待分类向量。默认第一行是横坐标，第二行是纵坐标。labelVector : np.matrix一个由+1和-1组成的1行m列矩矩阵，作为dataMatrix各列向量的标签。A : np.matrix一个1行m列矩阵，代表各待分向量对应的条件系数，非零系数对应支持向量。默认为空，即不画分界线。b : float一个1行1列矩阵，代表决策平面的截距。默认为空，即不绘制分界线。kernel : function核函数。用于在低维空间计算映射到高维空间的向量内积。为空则默认为原向量空间内积。sigma : float核函数的调控参数。越小映射的空间维度约高。Returns-------nullNotes-----只能绘制二维向量散点图。Examples--------X = np.matrix('1, 1; 4, 3; 3, 3').T Y = np.matrix('-1, 1, 1')A = np.matrix('0.25, 0, 0.25')b = np.matrix('-2')showDataSet(X, Y, A, b)'''n, m = np.shape(dataMatrix)# 初始化待分向量维度和向量个数if not n == 2: # 校验向量维度只能是2维raise Exception("only 2-dimension vectors can be darwn")plusColumNum = []# 正向量列号集合minusColumNum = [] # 负向量列号集合for i in range(m):if labelVector[0,i] > 0: # 如果标签为正plusColumNum.append(i) # 将列号归入正向量列集合else:minusColumNum.append(i) # 否则将列号归入负向量列集合plusData = dataMatrix[:, plusColumNum] # 由正向量列号获取组成正向量矩阵minusData = dataMatrix[:, minusColumNum] # 由负向量列号获取组成负向量矩阵plt.figure()plt.axis('scaled') # 横纵坐标尺度一致（即使在窗口缩放时）x_min = min(min(dataMatrix[0].tolist()[0])-1, 0) # X轴下限x_max = max(max(dataMatrix[0].tolist()[0])+1, 0) # X轴上限y_min = min(min(dataMatrix[1].tolist()[0])-1, 0) # Y轴下限y_max = max(max(dataMatrix[1].tolist()[0])+1, 0) # Y轴上限plt.xlim([x_min, x_max]) # 设置x轴坐标系范围plt.ylim([y_min, y_max]) # 设置y轴坐标系范围plt.scatter(plusData[0].tolist()[0], plusData[1].tolist()[0]) #正向量散点图plt.scatter(minusData[0].tolist()[0], minusData[1].tolist()[0]) #负向量散点图# 移动坐标系ax = plt.gca() # 获取当前坐标系ax.spines['right'].set_color('none') # 右边框设置成无颜色ax.spines['top'].set_color('none') # 上边框设置成无颜色 ax.spines['bottom'].set_position(('data',0)) # x轴在y轴，０的位置ax.spines['left'].set_position(('data',0)) # y轴在x轴，０的位置# 绘制分界线if b is not None: # 如果条件系数向量非空x = np.linspace(x_min, x_max, 100)y = np.linspace(y_min, y_max, 100)meshX,meshY = np.meshgrid(x,y) # 对向量显示区域网格化vecX = matlib.zeros((2, 1)) # 初始化决策面方程自变量向量Z = matlib.zeros((100, 100)) # 初始化曲面高度坐标for i, item_x in enumerate(x):for j, item_y in enumerate(y):vecX[0] = item_xvecX[1] = item_ymatX = vecX * matlib.ones((1, m))if callable(kernel): # 如果给定了可调用的核函数# 使用核函数计算每个网格点的坐标与所有乘过条件系数和标签值后的待# 分向量的高维内积，再加截距即得该网格点上的曲面高度Z[j,i] = kernel(dataMatrix, matX, sigma).T * np.multiply(labelVector, A).T +belse: # 否则直接在当前二维空间做内积# 此处务必注意i是行号，其实是y坐标；j是列号，其实是x坐标Z[j,i] = (dataMatrix.T * vecX).T * np.multiply(labelVector, A).T +bplt.contour(meshX, meshY, np.array(Z), [0], colors='r')# 标注支持向量和异常向量if A is not None and np.shape(A)[1] == m: # 如果条件系数向量非空svOder = []# 特殊向量序号集合for i in range(m):if abs(A[0,i]) > 0.01:# 条件系数显著大于0的才是特殊向量svOder.append(i)svSet = dataMatrix[:, svOder] # 特殊向量集合plt.scatter(svSet[0].tolist()[0], svSet[1].tolist()[0], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')

（三）线性分割效果

import numpy as npimport numpy.matlib as matlibimport numpy.random as randomimport svnSimpleimport showDataSet# 设定待分类点的坐标及其标签plusX = np.matrix(random.standard_normal((20, 2))).T + np.matrix('1; 1')plusY = matlib.ones((1,np.shape(plusX)[1]))minusX = np.matrix(random.standard_normal((25, 2))).T + np.matrix('5; 5')minusY = matlib.ones((1,np.shape(minusX)[1]))*-1X1 = np.c_[plusX, minusX]Y1 = np.c_[plusY, minusY]# 使用SVM模型给出决策平面的法向量和截距，以及条件系数向量A, b = svnSimple(X1, Y1, 50, 50)# 绘制待分向量散点图和分界线showDataSet(X1, Y1, A, b)

效果图如下：

（四）非线性分割效果

import numpy as npimport numpy.matlib as matlibimport numpy.random as randomimport svnSimpleimport showDataSet# 设定非线性待分类点的坐标及其标签（同心圆）plusX = np.matrix(random.standard_normal((40, 2))).TplusY = matlib.ones((1,np.shape(plusX)[1]))minusXright = np.matrix(random.standard_normal((10, 2))).T + np.matrix('7; 0')minusXleft = np.matrix(random.standard_normal((10, 2))).T + np.matrix('-7; 0')minusXup = np.matrix(random.standard_normal((10, 2))).T + np.matrix('0; 7')minusXdown = np.matrix(random.standard_normal((10, 2))).T + np.matrix('0; -7')minusXru = np.matrix(random.standard_normal((10, 2))).T + np.matrix('5; 5')minusXrd = np.matrix(random.standard_normal((10, 2))).T + np.matrix('5; -5')minusXlu = np.matrix(random.standard_normal((10, 2))).T + np.matrix('-5; 5')minusXld = np.matrix(random.standard_normal((10, 2))).T + np.matrix('-5; -5')minusX = np.c_[minusXright, minusXleft, minusXup, minusXdown, minusXru, minusXrd, minusXlu, minusXld]minusY = matlib.ones((1,np.shape(minusX)[1]))*-1X2 = np.c_[plusX, minusX]Y2 = np.c_[plusY, minusY]# 使用SVM模型和高斯核函数给出非线性分割线的条件向量和截距A, b = svnSimple(X2, Y2, 300, 100, gaussKernel, 3)# 绘制待分向量散点图和分界线showDataSet(X2, Y2, A, b, gaussKernel, 3)

效果图如下：