支持向量机（Support Vector Machine）

支持向量机

linear regression ， perceptron learning algorithm ， logistics regression都是分类器，咱们可使用这些分类器作线性和非线性的分类，好比下面的一个问题：

这里的每一条线都是能够把这个平面分开的，支持向量机要作的就是要在这些能够选择的直线中选择一条最好的直线来做为分类的直线。再给一个简单的解释，好比下面的三个图片，圆圈区域越大，说明这条直线对这些点放错的容忍度就越高：

##①超平面 介绍SVM以前，先来看超平面的概念： 其实超平面就是用于分割当前维度的一个空间。好比一维能够用一个点来进行分割，二维用一条线来进行分割，那么这些点和线就叫作“超”平面。加双引号是由于他们其实并非正在的超平面，由于超平面是要求大于三维的。 因此四维空间里面的超平面就是三维。好比在二维的空间里的超平面就是一条直线：

三维里面的超平面：

(其实这里的应该都不能叫超平面，由于超平面是三维以及三维以上的) 咱们把a ， b ， c看作是W0 , W1 , W2...，把x , y , z看作是x1 , x2 , x3，那么就有：

而W向量就是这个平面的法向量，咱们要求的就是法向量，证实一下： 以二维为例：

相减获得：

而[ (X_0 - X_1) , (Y_0 , Y_1)]这个点就在这个平面上，因此获得了：

因此W就是平面的法向量，这上面的X表明的是这个平面而不是点。

##②函数间隔的最大化 刚刚说到支持向量机也不是找超平面了，而是找最好的超平面，也就是对于点的犯错的容忍度越大越好，其实就是函数间隔越大越好：

右边的明显要好过左边的，由于左边的可犯错空间大啊。因此咱们要寻找的就是最大最肥的超平面——hperplane。

函数的间隔最大，**其实就是距离直线距离最短的点离直线的距离要最大。**因此先要知道直线的距离怎么求： 首先咱们假设X0在平面上的投影是X1，则确定有法向量W垂直于X0X1,：

又由于：

右由于X1在平面上的，前半部分就是-b了，能够写成：

和上面那条式子相等就获得：

这就是咱们要求的点到直线的距离了。而若是这个hperplane是正确的话，那么全部点的分类都是对的，那么咱们就默认他是对的，因而有：

这里能够相乘的条件是，咱们默认label正确的是1错误的是-1，若是你的错误是0正确是1的话公式是不一样的。乘上一个Y首先是能够去掉绝对值，使得函数变得可微，另外乘上以后函数值的绝对值也不会有变化，使得求解更加方便。 因此，最后的咱们的优化目标就是这样了：

里面的minimize是指找到距离hperplane最小距离的点，最外面就是挑选一个最好的W,b使得这个距离最小的点距离hperplane是最大的。

##③目标函数的化简 对于上面的式子，注意看到里面的那个式子：

举一个例子：

咱们代入(4,3)这个点，获得19，彷佛这个数字太大了，咱们不想要他这么大，咱们两边同时除去19，这个时候咱们的超平面就变成了：

数字是边了，可是这个超平面仍是在这个位置，因此能够认为超平面是没有变化的，这就证实了咱们上面那个式子：

是能够经过对w,b进行放缩而把左边的结果放大的无限多倍的。既然这样，那这个东西留着有什么意义，直接放缩到1就能够了，因而咱们把他放缩到1，也就是最小值是1。其实等于1都是差很少的，由于最小值以后都是1，因而就是有了：

那么target fomula就能够变成这样：

放缩以后并非就完了，这个是要加入当作条件来使用的。对于这个：

事实上咱们不太喜欢化简这样的，找最大化的不就是找最小化的W吗？找最小化的W不就是找最小化的W*W^T吗？不如再加个1/2? 因此问题就变成了：

为何要加上1/2呢？实际上是为了后面求导的时候方便化简的，可是这样对结果又没有什么影响。而W变成平方其实就是用上凸优化，由于平方以后就个凸函数了，这样变换一样对于最优结果没有任何影响。 因此最后要优化的结果：

##④Dual problem and KKT condiction 对于上述有条件的最优化问题，天然就要用上lagrange乘子法了。

右边的约束条件是要小于等于0的，α ≥ 0，只不过前面是符号因此转一下而已。 到这里，其实只要把等式扔到Quadratic Programming里面直接计算就行了。下面的步骤其实都是围绕解决这个优化问题展开的。

---

---

先在这停顿一下，咱们考虑一下：

####⑴SVM的机器学习可行性问题： 首先先来观察一下这个式子，感受似曾相识。他和L2 regularization很像，对于L2 regularization，首先是先要计算Ein的最小值，所谓的Ein就是该模型在当前的训练数据上犯错误的指望值。而后再正则化，因此L2是Minimizing Ein and Regularized L2 Paradigms；而支持向量机正好相反，他是先假设我这个平面是分类正确的，而后minimize W方：

后面咱们会用这个结论的。 回顾一下机器学习要解决的问题： ①Ein ≈ Eout Ein刚刚解释过了，Eout就是这个model在全局所犯的错误，Ein ≈ Eout就是要求这个model是能够反映全局的，若是不能反映，那就是过拟合了。 ②Ein ≈ 0 这个就是训练错误要接近于0，在这一步就容易发生过拟合的现象了。 而Ein ≈ Eout，也就是泛化能力是被VC dimension限制的，也就是说，越复杂的模型他的VC dimension越复杂。也就是VC bound右边的Ω会很大，VC bound就会很大，致使Ein 远远小于Eout了，由于复杂的模型意味着更加小的Ein。 再提一下，VC dimension就是break point - 1获得的。 若是是这样的话，那么正常来讲SVM的VC dimension也会很大啊，由于他的W是和数据维度相关的，数据多少维，那么W就多少个，而W表明的是自由度，一般也就表明这VC dimension， 可是SVM的效果仍是很好，为何呢？是什么东西限制着SVM的VC dimension？ 咱们来看一个例子： 在一个圆上，有三个点，你想找到一条能够分开的直线，能够获得VC dimension是3（以前有同窗看到在一个圆上分类他的VC dimension是无限的的，这是由于有无数多个点给你玩，这里就三个点，无限你又用不了，因此就只能是三个了啦），可是若是加上限制条件，这条线宽是5，那么VC dimension就是0了，由于没有地方塞进去。因此若是是large margin，VC dimension ≤ 3的。如图：

因此，large margin就是SVM的VCdimension的限制条件，致使它的分类效果很好，VC dimension小了天然泛化能力就行了，这里就解决了Ein ≈ Eout的问题，Ein ≈ 0这就是咱们后面要用凸优化来解决的问题了。

---

---

回到正题：

如何优化咱们的target function？ 上面的讨论咱们已经获得了VC dimension ≤ d + 1，**咱们悲观估计一下，就算SVM的VC dimension是d + 1了，加上d就是数据维度，加上的1是偏置值b。那么若是数据维度很大的话，计算复杂度是很高的，另外，如今咱们所研究的SVM仍是linear separable的，若是是来个nonlinear transform，数据维度就更加大了，再加上通常数据数量都是很的，时间会很长。**因此咱们要想一个方法来把数据维度d转移到其余地方去或者之间丢了。 而Daul problem刚好就能够解决。 回到origin target function： 咱们须要最小化：

对原函数作一些变换：

当这个点是违反了条件的时候，那么约束条件就会 > 0，α > 0，再maximum α那么就是无限大了，这个时候就不多是最小值，不可能取他，由于是无限大了。 当这个点是不违反的时候，那么约束条件就会 < 0，α > 0，再maximum α约束条件就是0，minimize w，b以后就仍是原来的target function。 因此变换以后：

变换以后的问题 == origin target function

---

---

再次停顿一下，考虑一下KKT条件是什么：

####⑵KKT 条件的引出 对于上述装换过的target function，有以下性质：

那么对于任何的条件都会有左边≥右边的。左边再加上一个w,b取最小：

一样是大于右边的全部状况，那么天然了，我右边再加上一个取α的最大值也是能够的：

而在右边的咱们把条件minimize和maximum调换过的式子就叫作Daul Problem。**因此，原问题是≥对偶问题的。**那么有没有什么办法能够把原问题的求解转换成对偶问题的求解呢？

---

---

####⑶KKT 条件的简单证实

对偶的意思就是存在一个最优的解使得两边的等式成立。因此咱们假设有一个W和B是最优的，那么有：

而最后能够看到求出来的解正是咱们要求的f(W)原目标函数，而原式子：

代进去也将是这个结果，由于maximum以后ag(x) = 0，因此本质上这个函数仍是求f(W)的最小值。因此对偶式子和原式在结果上是没有差异的。根据上面的式子，咱们本质就是要求

的最小值，固然这里的W,B要替换成原来的变量w，b了。求最小值天然就是求梯度为0了，因此w,b梯度为0的条件就有了。还有一个ag(x) = 0的条件，这个实际上是前提条件：

咱们以前说这个式子是等同目标函数的，既然要等同天然是要把后面的g(x)和h(x)消去啊！而h(x) = 0原本就消去了，而g(x) < 0，求最大必然就ag(x) = 0了，由于只有这个条件，才能消去后面的ag(x)把这个minimum maximum式子变成minimumf(w)的式子。因此再加上先前的几个拉格朗日条件就组成了KKT条件了。 因此KKT condition就是：

最后的几个条件实际上是lagrange乘子法的条件。

---

---

回到正题，既然咱们知道了能够利用KKT condition把origin target function转换到daul problem 来求解，那么上面这个问题咱们尝试用KKT条件求解一下： 首先对w，b求偏导：

把结果代回到dual problem：

因此最后咱们的target function就变成了这样。 最后咱们能够用QP对这个问题进行求解，求出了α以后，咱们随便取一个α是非0的，也就是＞0的解，这时候利用α*g(x) = 0的条件获得

b就求出来了，对于w直接代换上面的公式就行了。 而当α>0，由上面的公式能够获得这个点就刚恰好是在边界上，而这些点就叫作support vector，支持向量的点。咱们的拟合直线也将会由着些点肯定，其余不是support vector的点α就是0。

---

---

又停顿一下，咱们对这个式子思考一下：

####⑷为何咱们须要dual problem 其实这个最优问题用普通的QP求解也是能够的，可是若是咱们的数据维度很大，而通过feature transform以后维度就好更大了，这样就会致使VC dimension会增大，特别是后面用的多项式核，RBF核等等。通过对偶问题KKT条件的变换以后，咱们的目标式子：

转换成对偶问题以后，变量个数是N个，约束条件也是N个，于VC dimension就没有了关系，从某种意义上是简化了计算复杂度。其实计算复杂度仍是没有变，只是把维度的计算提高到了变量之间点的內积罢了。将原始SVM转化为对偶问题，本意是在非线性变化，进行特征转换后，若是d’很大，为了简化计算，消除d’的影响。进一步引入Kernel SVM，根本上解决上述问题。注意了，这里只是从某个角度看确实是消除了d维度的影响，实际上并无消失，只是转移到了计算內积里面而已。

---

---

回到正题，咱们的target function：

至于这个α怎么求，后面会用SMO算法求解。 到这里linear SVM就算结束了。

这就是分类函数。

---

---

再停顿一下，**什么是支持向量点，为何非支持向量的点α = 0？**这里仅仅思考linear SVM，若是是soft margin又不同了。

####⑸支持向量

若是是支持向量，他的function margin是1；而对于很多支持向量的点，function margin > 1，因此右边是负数，为了知足最大，因此α只能为0了，因此非支持向量的点α就是0。

---

---

##⑤kernel Support Vector Machine 回到正题，刚刚只是讲了linear SVM，是对于linear separable有效而已，若是是linear inseparable呢？好比一个圆形，这样就玩不了。记得以前linear regression和logistics regression讲到过一个feature transform，若是是非线性的咱们能够映射到其余维度进行解决，好比最多见的polynomial transform，可是这样问题来了，刚刚不是才把维度d转移到內积吗？（用dual problem的KKT condition）在来个feature transform那就是φ(x0)φ(x1)了，维度就更大了。

好比polynomial：

二项式的是这样的，注意到中间好像多了一个X1Xd，这是为了后面计算方便而已。两个作內积：

能够看到，最后的转换就只和原始空间有关系而已，对于转换只后的z空间的维度没有关系。好比x空间是2维的，为了解决nonlinear problem，咱们映射到了z空间，在z空间里面维度确定会比在x空间的原始维度要大，而最后用z空间作內积咱们就只须要拿x空间的原始维度就行了，由于咱们能够先內积再升维，而不是先升维再內积。

这种就叫作核函数了。 最后的分类函数用kernel function替代：

刚刚所讲的就是核函数的一种—— polynomial kernel function

加上几个参数，γ就是它的参数了，最后化简一下：

虽然都是二次转换，对应到同一个z空间。可是，若是他们的γ系数不一样，内积就会不同，那么就表明有不一样的距离，最终可能会获得不一样的SVM margin。因此，系数不一样，可能会获得不一样的hperplane。 看一下γ系数对于hperplane的影响：

使用高阶的polynomial kernel的话，获得的Support Vector数量不会太多，分类面不会太复杂，防止过拟合。一样也避开了对升维以后维度的依赖。

接下来介绍另一种kernel function——Gaussion kernel function 刚刚介绍的Q阶多项式是有限维度的，若是是无限维度的能不能经过kernel来简化计算？？有一个无限维的kernel function——Gaussion kernel

这和咱们以前见的有些不一样，只是去掉了下面的方差而已，方差是定值没有什么太大的影响。逆推看看它的维度是多少：

推出来后面的维度是无限个（中间用的是Taylor展开，由于e的特殊求导性质能够简化）。

分类函数就出来了。 可是核函数的过拟合仍是有一点严重的：

γ对于核函数的影响有点大。若是取值很大的话最后就会造成一个一个的小圈圈把那些点圈起来。

---

---

又得停顿一下，思考一下核函数的意义以及他们之间的对比：

####⑹Comparison of Kernels Polynomial Kernel的hyperplanes是由多项式曲线构成。优势：阶数能够灵活设置，更贴近实际分布；缺点：当Q很到的时候，若是kernel里面的值算出来是<1，那就基本接近0了，大于1就会变得很大，增长计算复杂度。并且参数过多，难以选择合适的值。

Gaussan Kernel的优势是边界更加复杂多样，能最准确地区分数据样本，数值计算K值波动较小，并且只有一个参数，容易选择；缺点是因为特征转换到无限维度中，w没有求解出来，计算速度要低于linear kernel，并且可能会发生过拟合。mysterious——no w；slower；too powerful。

以前说过经过对偶问题，咱们的把数据维度转移到了內积上，因此从某一方面来看咱们确实是作到了简化计算复杂度，可是实际上內积仍是属于一个很大的计算。 因此核函数的功能之一，就是简化计算，把升维和计算內积合在了一块儿，减小计算复杂度。把计算步骤结合在了一块儿，以前是先映射再计算內积，如今是一块儿作了。核函数的功能之二，就是能够很好的计算两个样本点的类似性，即內积。既然是表明类似性，咱们可不可使用其余的核函数呢？或者本身建立一个，好比欧氏距离，余弦距离等等？答案是不行。 先来看一下kernel的矩阵：

这有点像以前的协方差矩阵，只是没有减去均值，因此对称半正定是基本性质了。因此天然，咱们本身建立或选择的时候也要选择 ①symmetric对称②positive semi-definite 半正定。这也是核函数有效性的判断。

---

---

回到正题，刚刚只是讲了一下对核函数的理解。 ##⑥Soft-Margin Support Vector Machine 上面应用到的Gaussion Kernel貌似仍是会出现过拟合，并且仍是蛮严重的，这说明large margin已经限制不了Gaussion kernel了，咱们须要找其余方法来处理这个问题。 以前有一个比较简单的算法——perceptron learning algorithm 这个算法对于nonlinear problem有一个很好的处理方式，咱们不要求必定要分类正确，咱们只要求找到一个错误最少的分类就能够了。因此他的function是这样：

不正确的就加个1，最小为止。SVM也能够用这种方法来限制。

加上一个条件，C参数就是对于这些错误惩罚度是多少，条件也变了，正确的≥ 1，错误的无论他。无论是小错仍是大错。 整合一下：

这个式子其实没有什么用， 首先不是线性的，用不了二次规划，更不要说对偶这些了，其次大错小错都是同等对待，connot distinguish small error and large error。 对于上述方案继续修正： 咱们采用一个ξ做为一个犯错程度，程度越大，惩罚越大。惩罚就是这个式子数值会变大，而后SVM要花更多的力气去处理。

接下来就是对偶问题的推导，和以前的hard其实差很少的，lagrange 乘子法加对偶条件：

一样，KKT条件：

C - α = β 因此有：0 < α < C 其余的基本一致：w求导为0：

b求导：

接下来就是求b了：

求b的公式里面有一个矛盾的地方，就是咱们要求b首先得要求出来ξ的值，可是ξ的值也只有b的公式能够求的处理，因此这就有一个鸡生蛋蛋生鸡的问题。因此咱们口语去掉这个ξ。咱们刚刚用到的是拉格朗日乘子法，后面的β(-ξ)是一个仿射函数，仿射函数有β(-ξ) = 0的性质，因此把β代换一下就获得了上图的公式。那么去掉ξ就是等于0了，那么就只有C-α不等于0才有啊，因此当这个α ∈ （0 ， C）的时候就有ξ为0，然后面咱们会讲到当α∈（0，C）的时候这个点实际上是支持向量的点。这样就能够求出了b。 接下来看看C取值：

直接从我之前在CSDN里面写过的拷贝过来了。

接下来看一下一个比较重要的东西： physical significance of α

为何βξ = 0？缘由和前一个公式是同样的，由于要取最大值，因此这里要等于0，β ≥ 0，而实际公式是negative ξ，因此乘上去要是0才能有最大；第二，若是不是等于0就不等因而原问题的求解了，不等于0就无故端多了一个inequality，和原问题不对等了。以后才能进行daul problem的转换。 咱们主要是从上面这两个公式来看当α取值不一样的时候对应的物理意义。

当α = 0，得ξ = 0，这个点就是没有放错的点，由于ξ = 0，不须要容忍。而α = 0，因此不是支持向量机的点，因此表明的就是在bound外而且分类正确的点。

当α∈（0，C），仍是获得ξ = 0，这时候就不同了，尚未错误的点，可是第一条式子括号里面等于0了，意味着就是在bound上的点，那么就是支持向量点了。

当α = C，不能肯定ξ是否是0了，

，表示就是错了多少，这种有两种状况，一种是分类正确了，可是距离太近；或者是分类错了。 当α > C，不存在的，上面都限制了。

理一下整个思路。 ①找到最好的hperplane，最宽的那个。 ②获得target function ③发现feature transform以后维度对于计算机复杂度有很大影响，用dual problem转移到內积处理 ④转移以后发现仍是复杂度在的，引出了kernel function ⑤发现kernel function仍是有overfitting的状况，因而又引入了soft margin

---

---

在讲SMO算法以前，先讲一下对于error function的理解：

####⑺对于SVM error function的理解 咱们把SVM换一种形式。对于ξ，其实他是每个点距离边界有多远，一种是violating margin，即不知足y(wTz + b) ≥ 1，那么ξ就能够表示成：1 - y(wTz + b) > 0。第二种状况就是not violating margin，即这个点在边界以外，就是知足上述公式了，这个时候ξ就是0，咱们整合一下： ξ = max ( 1 - y(wTz + b) , 0 )，代换进原来的支持向量机公式：

这个就是支持向量机的error function，先预判了Ein = 0，也就是全对的状况，前面有说到。 这个function有点像咱们以前所学的L2 lost function：

这和logistics regression的L2范式的cost function很类似。

其实就差很少是同样的，没有什么差异，可是既然是相同的为何不用这种方法呢？两个缘由，一个是这种无条件的最优化问题没法经过QP解决，即对偶推导和kernel都没法使用；另外一个是这种形式中包含的max()项可能形成函数并非到处可导，这种状况难以用微分方法解决。

对比发现，L2 regularization和soft margin SVM形式是同样的，两个式子λ和C是互相对应的。soft marginSVM里面的large margin就对应着L2 regularization里面的short w，都是让hypothesis set能够简单点。λ和C也是互相对应，λ大，w就小，正则化的程度就越大；C小，Ein就大，响应这个margin也会打，因此增大C和减少λ是一个意思，因此large margin等同于regularization，都是防止过拟合做用的。

若是是按照咱们以前的err0/1，正确为1，错误就是0，那么有：

能够看到SVM他是大于err0/1的，根据VC bound理论是能够用来代替err0/1分类的。 后面再加上logic function的cost function：

而这个几乎就是和L2-regularized logistic regression同样的。Logistic Regression和Soft-Margin SVM都是在最佳化err0/1的上界而已。能够看出，求解regularized logistic regression的问题等同于求解soft-margin SVM的问题。

####⑻损失函数 常见的损失函数： err0/1：

此时soft margin就是这样了，大于0就是1小于就是0。 不敏感损失函数 —— hinge lost function

还有对数损失函数交叉熵等等。logistics用的是交叉熵，SVM就是用的hinge lost function。支持向量机就是一个结构风险最小化的近似实现，结构风险至关于指望风险(Eout)的一个上界，它是经验风险（Ein）和置信区间(Ω模型复杂度)的和，经验风险依赖于决策函数f的选取，可是置信区间是，F的VC维的增函数，二者是矛盾的。矛盾体如今：当VC维数变大的时候能够选到更好的f使得经验风险比较小，可是此时的置信区间比较大。这就是对应了VC bound理论。还好去听了台湾大学林轩宇老师课程，对这些机器学习理论基础有了解。

---

---

回到正题，开始SMO算法。 ##⑦SMO算法 target function：

刚刚咱们知道怎么求w,b，可是那是在知道了α的前提下，如今就来求α。 基本思路： 选择两个变量，固定其余变量，针对两个变量构建一个二次规划问题。每次针对两个变量来求解目标函数的最小值，求解完后，继续寻找新的变量求目标函数，在每次寻找新α的过程当中，目标函数将进一步获得优化，直到全部的αi更新完了。而对于α的选取，一个是违反KKT条件最严重的那一个，另外一个由约束条件自动肯定。

首先，假设咱们选取了两个变量α1，α2，固定其余变量以后：

因此只要求出α2，α1就知道了。

原目标函数化简以后：

K11指的就是x1和本身自己作核函数。因为咱们已经固定了除了α1和α2，因此天然其余的常量咱们能够去掉了，不如优化w+1，和优化w是同样的，去掉固定常数项就留下了上图的公式。

别忘了条件，条件是后面求解的关键。 首先咱们要获得α1，α2的范围。

由着两个约束条件限制。

因此有：

因此当咱们更新了α以后，咱们还要根据范围剪辑α才能够。

咱们假设：

剪辑范围：

再假设一个定值，也就是i = 3开始求和的：

目标式子：

用上面的vi代换以后：

求α2的话天然是求导了：

为0获得：

代入获得：

这里的化简有点麻烦：

手动证实一下。 用假设替换一下上面的式子：

就能够了。 SMO算法有两个要点：①α1的选择，违反KKT最严重的条件②α2的选择策略

---

---

很重要的问题，变量要怎么选择，后面会有例子证实。

####⑼变量的选择方式 SMO称选择第1个变量的过程为外层循环。外层循环在训练样本中选取违反KKT条件最严重的样本点，Violation of the most serious sample of KKT conditions。我第一次看这东西是懵逼的。可是仔细想一下，就是检测哪个样本是没有知足KKT的条件：

首先遍历全部0 < α < C的样本点，看看是否是知足的，若是没有载变量全部的。检测是否知足KKT。因此在SMO迭代的两个步骤中，只要α中有一个违背了KKT条件，这一轮迭代完成后，目标函数的值必然会增大。Generally speaking，KKT条件违背的程度越大，迭代后的优化效果越明显，增幅越大。 α1选完了天然就是选择第二个α了，第二个变量的选择叫作内存循环，咱们这里先用普通随机选择，看看效果如何。

---

---

##⑧算法实现——version 1 首先是导入各类各样的包和一个工具了：

import numpy as np  
import matplotlib.pyplot as plt  
import random  
import seaborn as sea  
import pandas as pd  


def get_positive_and_negative():  
  dataSet = pd.read_csv('Datas/LogiReg_data.txt', names=['V1', 'V2', 'Class'])  
  dataSet.Class[dataSet.Class == 0] = -1  
  dataSet = dataSet[60 : 80]  
  positive = dataSet[dataSet['Class'] == 1]  
  negative = dataSet[dataSet['Class'] == -1]  
  return positive , negative , dataSet  


def show_picture(positive , negative):  
  columns = ['V1', 'V2']  
  fig, ax = plt.subplots(figsize=(10, 5))  
  ax.scatter(positive[columns[0]], positive[columns[1]], s=30, c="b", marker="o", label="class 1")  
  ax.scatter(negative[columns[0]], negative[columns[1]], s=30, c="r", marker="x", label="class -1")  
  ax.legend()  
  ax.set_xlabel('V1')  
  ax.set_ylabel('V3')  
  plt.show()  

def load_data_set():  
  _ , _ , file = get_positive_and_negative()  
  orig_data = file.as_matrix()  
  cols = orig_data.shape[1]  
  data_mat = orig_data[ : , 0 : cols-1]  
  label_mat = orig_data[ : , cols-1 : cols]  
  return  data_mat , label_mat  

positive , negative , data = get_positive_and_negative()  
show_picture(positive , negative)  
print(data)  
复制代码

第一个是获得正负样本，而后显示，最后一个是加载数据，数据随便找一个就行了。

positive , negative , data = get_positive_and_negative()  
show_picture(positive , negative)  
复制代码

最后调用一些看看这些点是什么：

还有一些是对α的限制和一下工具函数：

''''' Generate a random number '''  
def select_jrand(i , m):  
  j = i  
  while(j == i):  
      j = int(random.uniform(0 , m))  
  return j  
  pass  

''''' restraint the α '''  
def clip_alpha(aj , H , L):  
  if aj > H:  
      aj = H  
  elif aj < L:  
      aj = L  
  return aj  
  pass  
复制代码

接下来就是实现支持向量机了：

def SVM(data_mat , class_label , C , tolar , max_iter):  

  data_mat = np.mat(data_mat)  
  label_mat = np.mat(class_label)  
  b = 0  
  m , n = np.shape(data_mat)  
  alphas = np.zeros((m , 1))  
  iter = 0  

  while iter < max_iter:  
      #做为迭代变化量 
      alpha_pairs_changed = 0  
      #做为第一个a 
      for i in range(m):  
          WT_i = np.dot(np.multiply(alphas , label_mat).T , data_mat)  
          f_xi = float(np.dot(WT_i , data_mat[i , :].T)) + b  
          Ei = f_xi - float(label_mat[i])  
          if ((label_mat[i]*Ei < -tolar) and (alphas[i] < C)) or ((label_mat[i]*Ei > tolar) and (alphas[i] > 0)):  
              j = Tools.select_jrand(i , m)  
              WT_j = np.dot(np.multiply(alphas , label_mat).T , data_mat)  
              f_xj  = float(np.dot(WT_j , data_mat[j , :].T)) + b  
              Ej = f_xj - float(label_mat[j])  
              alpha_iold = alphas[i].copy()  
              alpha_jold = alphas[j].copy()  

              if (label_mat[i] != label_mat[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 H == L:  
                  continue  

              eta = 2.0 * data_mat[i, :] * data_mat[j, :].T - data_mat[i, :] * data_mat[i, :].T - data_mat[j, :] * data_mat[j, :].T  
              if eta >= 0: continue  
              alphas[j] = (alphas[j] - label_mat[j]*(Ei - Ej))/eta  
              alphas[j] = Tools.clip_alpha(alphas[j], H, L)  
              if (abs(alphas[j] - alpha_jold) < 0.00001):  
                  continue  
              alphas[i] = alphas[i] + label_mat[j]*label_mat[i]*(alpha_jold - alphas[j])  


              b1 = b - Ei + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[i,:].T) +\  
              label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[i,:], data_mat[j,:].T)  
              b2 = b - Ej + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[j,:].T) +\  
              label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[j,:], data_mat[j,:].T)  
              if (0 < alphas[i]) and (C > alphas[i]):  
                  b = b1  
              elif (0 < alphas[j]) and (C > alphas[j]):  
                  b = b2  
              else:  
                  b = (b1 + b2)/2.0  
              print(b)  
              alpha_pairs_changed += 1  
              pass  
      if alpha_pairs_changed == 0:  
          iter += 1  
      else:  
          iter = 0  

  support_x = []  
  support_y = []  
  class1_x = []  
  class1_y = []  
  class01_x = []  
  class01_y = []  
  for i in range(m):  
      if alphas[i] > 0.0:  
          support_x.append(data_mat[i, 0])  
          support_y.append(data_mat[i, 1])  
  for i in range(m):  
      if label_mat[i] == 1:  
          class1_x.append(data_mat[i, 0])  
          class1_y.append(data_mat[i, 1])  
      else:  
          class01_x.append(data_mat[i, 0])  
          class01_y.append(data_mat[i, 1])  
  w_best = np.dot(np.multiply(alphas, label_mat).T, data_mat)  
  fig, ax = plt.subplots(figsize=(10, 5))  
  ax.scatter(support_x, support_y, s=100, c="y", marker="v", label="support_v")  
  ax.scatter(class1_x, class1_y, s=30, c="b", marker="o", label="class 1")  
  ax.scatter(class01_x, class01_y, s=30, c="r", marker="x", label="class -1")  
  lin_x = np.linspace(0, 100)  
  lin_y = (-float(b) - w_best[0, 0] * lin_x) / w_best[0, 1]  
  plt.plot(lin_x, lin_y, color="black")  
  ax.legend()  
  ax.set_xlabel("factor1")  
  ax.set_ylabel("factor2")  
  plt.show()  
  return b , alphas  
datamat , labelmat = dataSet.load_data_set()  
b, alphas = SVM(datamat , labelmat , 0.6 , 0.001 , 10)  
print(b , alphas)  
复制代码

首先传入的后面几个参数分别是惩罚力度，容忍度。比较重要的应该是这一句：

if ((label_mat[i]*Ei < -tolar) and (alphas[i] < C)) or ((label_mat[i]*Ei > tolar) and (alphas[i] > 0)):  
复制代码

这句话翻译过去就是yg(x) < 1 - ξ或者是y(g(x)) > 1+ξ。若是是小于，则这个点是离hperplane比较近，这时候这个点应该是等于C才对的；若是是大于了，也就是远大于边界了，那就是离边界很远了，可是α又大于0，离边界远意味着不是支持向量，因此α应该是0，因此能够改变。 后面的那些就是依据公式来的：

每一条都是对应公式写的。 最后就是打印了。

效果：

能够看到是极度不稳定。这是几个月前我实现的，后来如今我又从新实现了一个，用了一些改进方法。 为何会不稳定，我总结了几个缘由： ①没有缓存，更新慢，迭代次数不够 ②对于α2的选取没有很好的采起策略 ③对于n，也就是更新公式：

我没有判断是否是大于0的。n是什么东西呢？

他要是小于0意味着这个kernel matrix就不是半正定的了，K11 + K22 < 2K12；另外，这个n实际上是：

的二阶导数，小于0就不是凸函数了，哪来的凸优化。因此应该是更新的时候遇到这些状况致使不稳定的。 基于上面的缺点更换策略。 ##⑨算法实现——version 2 首先要改变的是加上一个缓存，用来保存Ei的值，使得计算更块。**其次就是α2的选择策略，在优化过程当中，会经过最大化步长的方式来得到第二个alpha值。第二步优化为，数据集全程扫描策略与在非边界alpha对中进行更新策略交替进行。**对于n，会进行判断是否是大于0，在这里是用-号的，因此n与咱们表达式上的是想反方向，因此是大于0。 首先仍是工具：

''' load data and define some tool function '''
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random

def loadDataSet(filename):
  ''' :param filename: :return dataset and label: '''

  dataset = []
  label = []
  fr = open(filename)
  for line in fr.readlines():
      lineArr = line.strip().split('\t')
      dataset.append( [np.float32(lineArr[0]) , np.float32(lineArr[1])] )
      label.append(np.float32(lineArr[2]))
  return dataset , label
  pass

''' select alpha2 randomly '''
def selectAlphaTwo(i , m):
  ''' :param i: :param m: :return: '''
  j = i
  while(j == i):
      j = int(random.uniform(0 , m))
  return j

def rangeSelectionForAlpha(aj , H , L):
  if aj > H:
      aj = H
  if L > aj:
      aj = L
  return aj
  pass

''' calculate Ei '''
def calculateEi(os , k):
  fxk = float(np.multiply(os.alphas, os.labels).T * (os.x * os.x[k, :].T)) + os.b
  Ek = fxk - float(os.labels[k])
  return Ek

''' put the Ei into the cache when calculate Ei '''
def selectj(i , os , Ei):
  maxk = -1
  maxDeltaE = 0
  Ej = 0
  os.eCache[i] = [1 , Ei]
  validEachlist = np.nonzero(os.eCache[: , 0].A)[0]
  if (len(validEachlist) > 1):
      for k in validEachlist:
          if k == i:
              continue
          Ek = calculateEi(os , k)
          deltaE = np.abs(Ei - Ek)
          if deltaE > maxDeltaE:
              maxk = k
              maxDeltaE = deltaE
              Ej = Ek
      return maxk , Ej
      pass
  else:
      j = selectAlphaTwo(i , os.m)
      Ej = calculateEi(os , j)
  return j , Ej
  pass

''' draw picture '''
def drawDataset(data , label , x = None , y = None , line = True , alphas = None , kernel = True):
  index_one = []
  index_negative_one = []
  for i in range(100):
      if label[i] == 1:
          index_one.append(data[i])
      else:
          index_negative_one.append(data[i])
  index_one = np.matrix(index_one)
  index_negative_one = np.matrix(index_negative_one)
  plt.scatter(index_one[ : , 0].tolist() , index_one[: , 1].tolist() , c = 'r' , marker='<' , label = 'class equal one')
  plt.scatter(index_negative_one[: , 0].tolist() , index_negative_one[: , 1].tolist() , c = 'b' , marker='x' , label = 'class equal negative one')
  if line == True:
      plt.plot(x , y)
      pass

  ''' draw the support vector,the point which the α not equal zero '''
  if line == True or kernel == True:
      a1 = []
      for i in range(len(alphas)):
          a = alphas[i]
          if a != 0:
             a1.append(data[i])
      a1 =  np.matrix(a1)
      print('The number of the support vector : ' , len(a1))
      plt.scatter(a1[: , 0].tolist(),a1[: , 1].tolist(), s=150, c='none', alpha=0.7,
                     linewidth=1.5, edgecolor='#AB3319' , label = 'support vector')

  plt.legend()
  plt.xlabel('X axis')
  plt.ylabel('Y axis')
  plt.show()

def updateEk(os,k):
  Ek = calculateEi(os,k)
  os.eCache[k]=[1,Ek]

if __name__ == '__main__':
  data , label = loadDataSet('../Data/testSetRBF.txt')
  drawDataset(data , label , line=False ,kernel=False)


复制代码

SMO算法惟一的一个类：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random
import KernelTransform
class optStruct:
  def __init__(self , dataMat , labels , C , toler):
      self.x = dataMat
      self.labels = labels
      self.C = C
      self.toler = toler
      self.m = np.shape(dataMat)[0]
      self.alphas = np.mat(np.zeros((self.m , 1)))
      self.b = 0
      self.eCache = np.mat(np.zeros((self.m , 2)))
      self.K = np.mat(np.zeros((self.m , self.m)))
      for i in range(self.m):
          self.K[: , i] = KernelTransform.kernelTrans(self.x , self.x[i , :] , kTup=('rbf' , 1.2))
      pass

if __name__ == '__main__':
  os = optStruct([1,2] , [3,4] , 1,1)
  a = os.alphas.tolist()[0][0] -  os.alphas.tolist()[1][0]
  print(max(1.0 , a))


复制代码

须要解释的应该只有selectj()了，这个是经过计算最大不长来选择α2的。 首先咱们假设最大不长是-1，由于相减有绝对值不多是negative；os.eCache是咱们的缓存的Ei，先把Ei存进去，1,表示这个数字不是0，这一步就是获得这个缓存里面全部有效（不为0）的Ei。判断获得的列表是否是有东西，没有就随机选择了。仍是再解释一下为何要这个创建表格吧！ 咱们在选择第一个α1的时候，选择的是在边界外的点，也就是非边界的点。 优先选择遍历非边界数据样本，由于非边界数据样本更有可能须要调整，边界数据样本经常不能获得进一步调整而留在边界上。因为大部分数据样本都很明显不多是支持向量，所以对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本获得调整时，遍历全部数据样本，以检验是否整个集合都知足KKT条件。若是整个集合的检验中又有数据样本被进一步进化，则有必要再遍历非边界数据样本。这样，不停地在遍历全部数据样本和遍历非边界数据样本之间切换，直到整个样本集合都知足KKT条件为止。以上用KKT条件对数据样本所作的检验都以达到必定精度ε就能够中止为条件。若是要求十分精确的输出算法，则每每不能很快收敛。因此在echa中缓存的第一次选出的α，由于咱们选出来的就是非边界上的点，α2选择的时候继续在上面遍历，虽然缓存是存了Ei，可是这个Ei不能直接用，由于那个是旧的值。因此α的迭代策略就是非边界和全局选取两种交替进行了。

以后就是正式的算法了：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random
import Tool
import smo_class
import KernelTransform
def innerL(i ,os):
  Ei = Tool.calculateEi(os , i)
  if ((os.labels[i]*Ei < -os.toler) and
      (os.alphas[i] < os.C)) or ((os.labels[i]*Ei > os.toler) and
                                 (os.alphas[i] > 0)):
      j , Ej = Tool.selectj(i , os , Ei)
      alphaIold = os.alphas[i].copy()
      alphaJold = os.alphas[j].copy()
      if (os.labels[i] != os.labels[j]):
          L = max(0 , os.alphas[j] - os.alphas[i])
          H = min(os.C , os.C + np.array(os.alphas)[j] - np.array(os.alphas)[i])
      else:
          L = max(0 , os.alphas[j] + os.alphas[i] - os.C)
          H = min(os.C , np.array(os.alphas)[j] + np.array(os.alphas)[i])
      if L == H:
          return 0
      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('η> 0，the kernel matrix is not semi-positive definite')
          return 0
      os.alphas[j] -= os.labels[j]*(Ei - Ej)/eta
      os.alphas[j] = Tool.rangeSelectionForAlpha(os.alphas[j] , H , L)
      Tool.updateEk(os , j)

      if (abs(os.alphas[j] - alphaJold) < 0.00001):
          print("j not moving enough")
          return 0
      os.alphas[i] += os.labels[j] * os.labels[i] * (alphaJold - os.alphas[j])
      Tool.updateEk(os , i)
      b1 = os.b - Ei - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.x[i, :] * os.x[i, :].T - os.labels[j] * \
           (os.alphas[j] - alphaJold) * os.x[i, :] * os.x[j, :].T
      b2 = os.b - Ej - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.x[i, :] * os.x[j, :].T - os.labels[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
      return 1
  else:
      return 0

def smo(data,labels,C = 0.6,toler = 0.001,maxIter = 40 , kernel = True):
  oS = smo_class.optStruct(np.mat(data),np.mat(labels).transpose(),C,toler)
  iter =0
  entireSet  = True
  alphaPairsChanged = 0
  while(iter < maxIter) and ((alphaPairsChanged >0) or (entireSet)):
      alphaPairsChanged = 0
      if entireSet:
          for i in range(oS.m):
              if kernel == True:
                  alphaPairsChanged += KernelTransform.innerL(i,oS)
              else:
                  alphaPairsChanged += innerL(i, oS)
          print("fullSet,iter: %d i: %d,pairs changed %d" %\
              (iter,i,alphaPairsChanged))
          iter +=1
      else:
          # 两个元素乘积非零，每两个元素作乘法[0,1,1,0,0]*[1,1,0,1,0]=[0,1,0,0,0]
          nonBoundIs = np.nonzero((oS.alphas.A > 0)*(oS.alphas.A < C))[0]
          for i in nonBoundIs:
              alphaPairsChanged += innerL(i,oS)
              print("nou-bound,iter: %d i:%d,pairs changed %d" % (iter,i,alphaPairsChanged))
          iter +=1
      # entireSet 控制交替的策略选择
      if entireSet:
          entireSet = False
      # 必须有alpha对进行更新
      elif(alphaPairsChanged == 0):
          entireSet = True
      print("iteration number：%d" % iter)
  return oS.b,oS.alphas
复制代码

entireSet就是交换策略的标志。貌似没有什么好说的。 以后就是执行函数这些了：

import Tool
import SMO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import KernelTransform
''' calculate w and draw the picture, the variable which the α not equal zero , we call support vector '''
def calculateW(alphas , data , labels):
  x = np.mat(data)
  label = np.mat(labels).transpose()
  m , n = np.shape(x)
  w = np.zeros((n , 1))
  for i in range(m):
      w += np.multiply(alphas[i] * label[i] , x[i , :].T)
  return w
  pass

if __name__ == '__main__':
  data, label = Tool.loadDataSet('../Data/testSet.txt')
  b,alphas = SMO.smo(data , label , kernel=False)
  w = calculateW(alphas , data , label)
  x = np.arange(0 , 11)
  print(w)
  y = (-b - w[0]*x)/w[1]
  Tool.drawDataset(data , label , x , y.tolist()[0] , line=True , alphas=alphas)

  data, label = Tool.loadDataSet('../Data/testSetRBF.txt')
  b, alphas = SMO.smo(data, label,kernel=True ,maxIter=100)
  svInd = np.nonzero(alphas.A > 0)[0]
  Tool.drawDataset(data, label,  line=False, alphas=alphas)







复制代码

有一个是kernel function的，先不用管。 效果：

圈起来的是支持向量点，好不少了。 ##⑩算法实现——version 3 kernel function加上，先看看原来的数据：

须要改的其实就是內积就能够了，处处看看哪里有內积就改改他，修改事后的innel和smo：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import Tool
def kernelTrans(X,A,kTup):
  m,n = np.shape(X)
  K = np.mat(np.zeros((m,1)))
  if kTup[0]=='lin':
      K = X*A.T
  elif kTup[0] =='rbf':
      for j in range(m):
          deltRow = X[j,:]-A
          K[j] = deltRow*deltRow.T
      K = np.exp(K/(-1*kTup[1]**2))
  return K

''' update the innel function '''
def innerL(i ,os):
  Ei = calculateEi(os , i)
  if ((os.labels[i]*Ei < -os.toler) and
      (os.alphas[i] < os.C)) or ((os.labels[i]*Ei > os.toler) and
                                 (os.alphas[i] > 0)):
      j , Ej = Tool.selectj(i , os , Ei)
      alphaIold = os.alphas[i].copy()
      alphaJold = os.alphas[j].copy()
      if (os.labels[i] != os.labels[j]):
          L = max(0 , os.alphas[j] - os.alphas[i])
          H = min(os.C , os.C + np.array(os.alphas)[j] - np.array(os.alphas)[i])
      else:
          L = max(0 , os.alphas[j] + os.alphas[i] - os.C)
          H = min(os.C , np.array(os.alphas)[j] + np.array(os.alphas)[i])
      if L == H:
          return 0
      eta = 2.0 * os.K[i, j] - os.K[i, i] - os.K[j, j]
      if eta >= 0:
          print('η> 0，the kernel matrix is not semi-positive definite')
          return 0
      os.alphas[j] -= os.labels[j]*(Ei - Ej)/eta
      os.alphas[j] = Tool.rangeSelectionForAlpha(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.labels[j] * os.labels[i] * (alphaJold - os.alphas[j])
      updateEk(os , i)
      b1 = os.b - Ei - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.K[i , i] - os.labels[j] * \
           (os.alphas[j] - alphaJold) *  os.K[i , j]
      b2 = os.b - Ej - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.K[i , j] - os.labels[j] * \
           (os.alphas[j] - alphaJold) * os.K[j , j]
      if (0 < os.alphas[i]) and (os.C > os.alphas[i]):
          os.b = b1
      elif (0 < os.alphas[j]) and (os.C > os.alphas[j]):
          os.b = b2
      else:
          os.b = (b1 + b2) / 2.0
      return 1
  else:
      return 0

''' updata the Ei '''
def calculateEi(os , k):
  fxk = float(np.multiply(os.alphas, os.labels).T * os.K[:, k] + os.b)
  Ek = fxk - float(os.labels[k])
  return Ek
def updateEk(os,k):
  Ek = calculateEi(os,k)
  os.eCache[k]=[1,Ek]
复制代码

刚刚那个执行函数其实已经包括了kernel的，因此直接就能够看到效果了：

用的是Gaussion kernel，不知道怎么作拟合，就把支持向量点圈出来就行了。 最后附上全部代码GitHub： github.com/GreenArrow2…