支持向量机（Support Vector Machine）支持向量机
支持向量机（Support Vector Machine）支持向量机
支持向量机
linear regression ， perceptron learning algorithm ， logistics regression都是分类器，我们可以使用这些分类器做线性和非线性的分类，比如下面的一个问题：
GV0SHYC3S{P{Q4QVB66UN6T.png
这里的每一条线都是可以把这个平面分开的，支持向量机要做的就是要在这些可以选择的直线中选择一条最好的直线来作为分类的直线。再给一个简单的解释，比如下面的三个图片，圆圈区域越大，说明这条直线对这些点放错的容忍度就越高：
0T0FN3C5MKTQE`DO{N}QX(C.png
①超平面
介绍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： https://github.com/GreenArrow2017/MachineLearning