支持向量机之SMO

[原文]http://blog.csdn.net/zouxy09/article/details/17291543

写的太吊了，不得不转，把我前大半个月看的论文都总结了，而且讲的非常好，最重要是有代码实现！

支持向量机

[e-mail]

http://blog.csdn.net/zouxy09

       机器学习算法与Python实践这个系列主要是参考《机器学习实战》这本书。因为自己想学习Python，然后也想对一些机器学习算法加深下了解，所以就想通过Python来实现几个比较常用的机器学习算法。恰好遇见这本同样定位的书籍，所以就参考这本书的过程来学习了。

       在这一节我们主要是对支持向量机进行系统的回顾，以及通过Python来实现。由于内容很多，所以这里分成三篇博文。第一篇讲SVM初级，第二篇讲进阶，主要是把SVM整条知识链理直，第三篇介绍Python的实现。SVM有很多介绍的非常好的博文，具体可以参考本文列出的参考文献和推荐阅读资料。在本文中，定位在于把集大成于一身的SVM的整体知识链理直，所以不会涉及细节的推导。网上的解说的很好的推导和书籍很多，大家可以进一步参考。

 

目录

一、引入

二、线性可分SVM与硬间隔最大化

三、Dual优化问题

       3.1、对偶问题

       3.2、SVM优化的对偶问题

四、松弛向量与软间隔最大化

五、核函数

六、多类分类之SVM

       6.1、“一对多”的方法

       6.2、“一对一”的方法

七、KKT条件分析

八、SVM的实现之SMO算法

       8.1、坐标下降算法

       8.2、SMO算法原理

       8.3、SMO算法的Python实现

九、参考文献与推荐阅读

 

一、引入

       支持向量机（SupportVector Machines），这个名字可是响当当的，在机器学习或者模式识别领域可是无人不知，无人不晓啊。八九十年代的时候，和神经网络一决雌雄，独领风骚，并吸引了大批为之狂热和追随的粉丝。虽然几十年过去了，但风采不减当年，在模式识别领域依然占据着大遍江山。王位稳固了几十年。当然了，它也繁衍了很多子子孙孙，出现了很多基因改良的版本，也发展了不少裙带关系。但其中的睿智依然被世人称道，并将千秋万代！

        好了，买了那么久广告，不知道是不是高估了。我们还是脚踏实地，来看看传说的SVM是个什么东西吧。我们知道，分类的目的是学会一个分类函数或分类模型（或者叫做分类器），该模型能把数据库中的数据项映射到给定类别中的某一个，从而可以用于预测未知类别。对于用于分类的支持向量机，它是个二分类的分类模型。也就是说，给定一个包含正例和反例（正样本点和负样本点）的样本集合，支持向量机的目的是寻找一个超平面来对样本进行分割，把样本中的正例和反例用超平面分开，但是不是简单地分看，其原则是使正例和反例之间的间隔最大。学习的目标是在特征空间中找到一个分类超平面wx+b=0，分类面由法向量w和截距b决定。分类超平面将特征空间划分两部分，一部分是正类，一部分是负类。法向量指向的一侧是正类，另一侧为负类。

        用一个二维空间里仅有两类样本的分类问题来举个小例子。假设我们给定了下图左图所示的两类点Class1和Class2（也就是正样本集和负样本集）。我们的任务是要找到一个线，把他们划分开。你会告诉我，那简单，挥笔一画，洋洋洒洒五颜六色的线就出来了，然后很得意的和我说，看看吧，下面右图，都是你要的答案，如果你还想要，我还可以给你画出无数条。对，没错，的确可以画出无数条。那哪条最好呢？你会问我，怎么样衡量“好”？假设Class1和Class2分别是两条村子的人，他们因为两条村子之间的地盘分割的事闹僵了，叫你去说个理，到底怎么划分才是最公平的。这里的“好”，可以理解为对Class1和Class2都是公平的。然后你二话不说，指着黑色那条线，说“就它了！正常人都知道！在两条村子最中间画条线很明显对他们就是公平的，谁也别想多，谁也没拿少”。这个例子可能不太恰当，但道理还是一样的。对于分类来说，我们需要确定一个分类的线，如果新的一个样本到来，如果落在线的左边，那么这个样本就归为class1类，如果落在线的右边，就归为class2这一类。那哪条线才是最好的呢？我们仍然认为是中间的那条，因为这样，对新的样本的划分结果我们才认为最可信，那这里的“好”就是可信了。另外，在二维空间，分类的就是线，如果是三维的，分类的就是面了，更高维，也有个霸气的名字叫超平面。因为它霸气，所以一般将任何维的分类边界都统称为超平面。

       好了。对于人来说，我们可以轻易的找到这条线或者超平面（当然了，那是因为你可以看到样本具体的分布是怎样的，如果样本的维度大于三维的话，我们就没办法把这些样本像上面的图一样画出来了，这时候就看不到了，这时候靠人的双眼也无能为力了。“如果我能看得见，生命也许完全不同，可能我想要的，我喜欢的我爱的，都不一样……”），但计算机怎么知道怎么找到这条线呢？我们怎么把我们的找这条线的方法告诉他，让他按照我们的方法来找到这条线呢？呃，我们要建模！！！把我们的意识“强加”给计算机的某个数学模型，让他去求解这个模型，得到某个解，这个解就是我们的这条线，那这样目的就达到了。那下面就得开始建模之旅了。

 

二、线性可分SVM与硬间隔最大化

      其实上面这种分类思想就是SVM的思想。可以表达为：SVM试图寻找一个超平面来对样本进行分割，把样本中的正例和反例用超平面分开，但是不是很敷衍地简单的分开，而是尽最大的努力使正例和反例之间的间隔margin最大。这样它的分类结果才更加可信，而且对于未知的新样本才有很好的分类预测能力（机器学习美其名曰泛化能力）。

      我们的目标是寻找一个超平面，使得离超平面比较近的点能有更大的间距。也就是我们不考虑所有的点都必须远离超平面，我们关心求得的超平面能够让所有点中离它最近的点具有最大间距。

       我们先用数学公式来描述下。假设我们有N个训练样本{(x₁, y₁),(x₂, y₂), …, (x_N, y_N)}，x是d维向量，而y_i∊{+1, -1}是样本的标签，分别代表两个不同的类。这里我们需要用这些样本去训练学习一个线性分类器（超平面）：f(x)=sgn(w^Tx + b)，也就是w^Tx + b大于0的时候，输出+1，小于0的时候，输出-1。sgn()表示取符号。而g(x) =w^Tx + b=0就是我们要寻找的分类超平面，如上图所示。刚才说我们要怎么做了？我们需要这个超平面最大的分隔这两类。也就是这个分类面到这两个类的最近的那个样本的距离相同，而且最大。为了更好的说明，我们在上图中找到两个和这个超平面平行和距离相等的超平面：H₁: y = w^Tx + b=+1 和 H₂: y = w^Tx + b=-1。

       好了，这时候我们就需要两个条件：（1）没有任何样本在这两个平面之间；（2）这两个平面的距离需要最大。（对任何的H₁和H₂，我们都可以归一化系数向量w，这样就可以得到H₁和H₂表达式的右边分别是+1和-1了）。先来看条件（2）。我们需要最大化这个距离，所以就存在一些样本处于这两条线上，他们叫支持向量（后面会说到他们的重要性）。那么它的距离是什么呢？我们初中就学过，两条平行线的距离的求法，例如ax+by=c1和ax+by=c2，那他们的距离是|c2-c1|/sqrt(x²+y²)（sqrt()表示开根号）。注意的是，这里的x和y都表示二维坐标。而用w来表示就是H₁:w₁x₁+w₂x₂=+1和H₂:w₁x₁+w₂x₂=-1，那H₁和H₂的距离就是|1+1|/ sqrt(w₁²+w₁²)=2/||w||。也就是w的模的倒数的两倍。也就是说，我们需要最大化margin=2/||w||，为了最大化这个距离，我们应该最小化||w||，看起来好简单哦。同时我们还需要满足条件（2），也就是同时要满足没有数据点分布在H₁和H₂之间：

      也就是，对于任何一个正样本y_i=+1，它都要处于H₁的右边，也就是要保证：y= w^Tx+ b>=+1。对于任何一个负样本y_i=-1，它都要处于H₂的左边，也就是要保证：y = w^Tx + b<=-1。这两个约束，其实可以合并成同一个式子：y_i (w^Tx_i + b)>=1。

      所以我们的问题就变成：

       这是个凸二次规划问题。什么叫凸？凸集是指有这么一个点的集合，其中任取两个点连一条直线，这条线上的点仍然在这个集合内部，因此说“凸”是很形象的。例如下图，对于凸函数（在数学表示上，满足约束条件是仿射函数，也就是线性的Ax+b的形式）来说，局部最优就是全局最优，但对非凸函数来说就不是了。二次表示目标函数是自变量的二次函数。

      好了，既然是凸二次规划问题，就可以通过一些现成的 QP (Quadratic Programming) 的优化工具来得到最优解。所以，我们的问题到此为止就算全部解决了。虽然这个问题确实是一个标准的 QP 问题，但是它也有它的特殊结构，通过 Lagrange Duality 变换到对偶变量 (dual variable) 的优化问题之后，可以找到一种更加有效的方法来进行求解，而且通常情况下这种方法比直接使用通用的 QP 优化包进行优化要高效得多。也就说，除了用解决QP问题的常规方法之外，还可以应用拉格朗日对偶性，通过求解对偶问题得到最优解，这就是线性可分条件下支持向量机的对偶算法，这样做的优点在于：一是对偶问题往往更容易求解；二者可以自然的引入核函数，进而推广到非线性分类问题。那什么是对偶问题？

 

三、Dual优化问题

3.1、对偶问题

       在约束最优化问题中，常常利用拉格朗日对偶性将原始问题转换为对偶问题，通过求解对偶问题而得到原始问题的解。至于这其中的原理和推导参考文献[3]讲得非常好。大家可以参考下。这里只将对偶问题是怎么操作的。假设我们的优化问题是：

min f(x)

s.t. h_i(x) = 0, i=1, 2, …,n

       这是个带等式约束的优化问题。我们引入拉格朗日乘子，得到拉格朗日函数为：

L(x, α)=f(x)+α₁h₁(x)+ α₂h₂(x)+…+α_nh_n(x)

       然后我们将拉格朗日函数对x求极值，也就是对x求导，导数为0，就可以得到α关于x的函数，然后再代入拉格朗日函数就变成：

max W(α) = L(x(α), α)

       这时候，带等式约束的优化问题就变成只有一个变量α（多个约束条件就是向量）的优化问题，这时候的求解就很简单了。同样是求导另其等于0，解出α即可。需要注意的是，我们把原始的问题叫做primal problem，转换后的形式叫做dual problem。需要注意的是，原始问题是最小化，转化为对偶问题后就变成了求最大值了。对于不等式约束，其实是同样的操作。简单地来说，通过给每一个约束条件加上一个 Lagrange multiplier（拉格朗日乘子），我们可以将约束条件融和到目标函数里去，这样求解优化问题就会更加容易。（这里其实涉及到很多蛮有趣的东西的，大家可以参考更多的博文）

 

3.2、SVM优化的对偶问题

       对于SVM，前面提到，其primal problem是以下形式：

       同样的方法引入拉格朗日乘子，我们就可以得到以下拉格朗日函数：

        然后对L(w, b, α)分别求w和b的极值。也就是L(w, b,α)对w和b的梯度为0：∂L/∂w=0和∂L/∂b=0，还需要满足α>=0。求解这里导数为0的式子可以得到：

       然后再代入拉格朗日函数后，就变成：

        这个就是dual problem（如果我们知道α，我们就知道了w。反过来，如果我们知道w，也可以知道α）。这时候我们就变成了求对α的极大，即是关于对偶变量α的优化问题（没有了变量w，b，只有α）。当求解得到最优的α*后，就可以同样代入到上面的公式，导出w*和b*了，最终得出分离超平面和分类决策函数。也就是训练好了SVM。那来一个新的样本x后，就可以这样分类了：

       在这里，其实很多的α_i都是0，也就是说w只是一些少量样本的线性加权值。这种“稀疏”的表示实际上看成是KNN的数据压缩的版本。也就是说，以后新来的要分类的样本首先根据w和b做一次线性运算，然后看求的结果是大于0还是小于0来判断正例还是负例。现在有了α_i，我们不需要求出w，只需将新来的样本和训练数据中的所有样本做内积和即可。那有人会说，与前面所有的样本都做运算是不是太耗时了？其实不然，我们从KKT条件中得到，只有支持向量的α_i不为0，其他情况α_i都是0。因此，我们只需求新来的样本和支持向量的内积，然后运算即可。这种写法为下面要提到的核函数（kernel）做了很好的铺垫。如下图所示：

 

四、松弛向量与软间隔最大化

       我们之前讨论的情况都是建立在样本的分布比较优雅和线性可分的假设上，在这种情况下可以找到近乎完美的超平面对两类样本进行分离。但如果遇到下面这两种情况呢？左图，负类的一个样本点A不太合群，跑到正类这边了，这时候如果按上面的确定分类面的方法，那么就会得到左图中红色这条分类边界，嗯，看起来不太爽，好像全世界都在将就A一样。还有就是遇到右图的这种情况。正类的一个点和负类的一个点都跑到了别人家门口，这时候就找不到一条直线来将他们分开了，那这时候怎么办呢？我们真的要对这些零丁的不太听话的离群点屈服和将就吗？就因为他们的不完美改变我们原来完美的分界面会不会得不偿失呢？但又不得不考虑他们，那怎样才能折中呢？

       对于上面说的这种偏离正常位置很远的数据点，我们称之为 outlier，它有可能是采集训练样本的时候的噪声，也有可能是某个标数据的大叔打瞌睡标错了，把正样本标成负样本了。那一般来说，如果我们直接忽略它，原来的分隔超平面还是挺好的，但是由于这个 outlier 的出现，导致分隔超平面不得不被挤歪了，同时 margin 也相应变小了。当然，更严重的情况是，如果出现右图的这种outlier，我们将无法构造出能将数据线性分开的超平面来。

      为了处理这种情况，我们允许数据点在一定程度上偏离超平面。也就是允许一些点跑到H₁和H₂之间，也就是他们到分类面的间隔会小于1。如下图：

       具体来说，原来的约束条件就变为：

        这时候，我们在目标函数里面增加一个惩罚项，新的模型就变成（也称软间隔）：

       引入非负参数ξ_i后（称为松弛变量），就允许某些样本点的函数间隔小于1，即在最大间隔区间里面，或者函数间隔是负数，即样本点在对方的区域中。而放松限制条件后，我们需要重新调整目标函数，以对离群点进行处罚，目标函数后面加上的第二项就表示离群点越多，目标函数值越大，而我们要求的是尽可能小的目标函数值。这里的C是离群点的权重，C越大表明离群点对目标函数影响越大，也就是越不希望看到离群点。这时候，间隔也会很小。我们看到，目标函数控制了离群点的数目和程度，使大部分样本点仍然遵守限制条件。

        这时候，经过同样的推导过程，我们的对偶优化问题变成：

       此时，我们发现没有了参数ξ_i，与之前模型唯一不同在于α_i又多了α_i<=C的限制条件。需要提醒的是，b的求值公式也发生了改变，改变结果在SMO算法里面介绍。

五、核函数

       如果我们的正常的样本分布如下图左边所示，之所以说是正常的指的是，不是上面说的那样由于某些顽固的离群点导致的线性不可分。它是真的线性不可分。样本本身的分布就是这样的，如果也像样本那样，通过松弛变量硬拉一条线性分类边界出来，很明显这条分类面会非常糟糕。那怎么办呢？SVM对线性可分数据有效，对不可分的有何应对良策呢？是核方法（kernel trick）大展身手的时候了。 

       如上图右，如果我们可以把我们的原始样本点通过一个变换，变换到另一个特征空间，在这个特征空间上是线性可分的，那么上面的SVM就可以轻易工作了。也就是说，对于不可分的数据，现在我们要做两个工作：

1）首先使用一个非线性映射Φ(x)将全部原始数据x变换到另一个特征空间，在这个空间中，样本变得线性可分了；

2）然后在特征空间中使用SVM进行学习分类。

       好了，第二个工作没什么好说的，和前面的一样。那第一个粗重活由谁来做呢？我们怎么知道哪个变换才可以将我们的数据映射为线性可分呢？数据维度那么大，我们又看不到。另外，这个变换会不会使第二步的优化变得复杂，计算量更大呢？对于第一个问题，有个著名的cover定理：将复杂的模式分类问题非线性地投射到高维空间将比投射到低维空间更可能是线性可分的。OK，那容易了，我们就要找到一个所有样本映射到更高维的空间的映射。对不起，其实要找到这个映射函数很难。但是，支持向量机并没有直接寻找和计算这种复杂的非线性变换，而是很智慧的通过了一种巧妙的迂回方法来间接实现这种变换。它就是核函数，不仅具备这种超能力，同时又不会增加太多计算量的两全其美的方法。我们可以回头看看上面SVM的优化问题：

      可以看到，对样本x的利用，只是计算第i和第j两个样本的内积就可以了。

       对于分类决策函数，也是计算两个样本的内积。也就是说，训练SVM和使用SVM都用到了样本间的内积，而且只用到内积。那如果我们可以找到一种方法来计算两个样本映射到高维空间后的内积的值就可以了。核函数就是完成这伟大的使命的：

K(x_i, x_j)=Φ(x_i)^T Φ(x_j)

       也就是两个样本x_i和x_j对应的高维空间的内积Φ(x_i)^T Φ(x_j)通过一个核函数K(x_i, x_j)计算得到。而不用知道这个变换Φ(x)是何许人也。而且这个核函数计算很简单，常用的一般是径向基RBF函数：

     这时候，我们的优化的对偶问题就变成了：

      和之前的优化问题唯一的不同只是样本的内积需要用核函数替代而已。优化过程没有任何差别。而决策函数变成了：

        也就是新来的样本x和我们的所有训练样本计算核函数即可。需要注意的是，因为大部分样本的拉格朗日因子α_i都是0，所以其实我们只需要计算少量的训练样本和新来的样本的核函数，然后求和取符号即可完成对新来样本x的分类了。支持向量机的决策过程也可以看做一种相似性比较的过程。首先，输入样本与一系列模板样本进行相似性比较，模板样本就是训练过程决定的支持向量，而采用的相似性度量就是核函数。样本与各支持向量比较后的得分进行加权后求和，权值就是训练时得到的各支持向量的系数α_i和类别标号的成绩。最后根据加权求和值大小来进行决策。而采用不同的核函数，就相当于采用不同的相似度的衡量方法。

       从计算的角度，不管Φ(x)变换的空间维度有多高，甚至是无限维（函数就是无限维的），这个空间的线性支持向量机的求解都可以在原空间通过核函数进行，这样就可以避免了高维空间里的计算，而计算核函数的复杂度和计算原始样本内积的复杂度没有实质性的增加。

       到这里，忍不住要感叹几声。为什么“碰巧”SVM里需要计算的地方数据向量总是以内积的形式出现？为什么“碰巧”存在能简化映射空间中的内积运算的核函数？为什么“碰巧”大部分的样本对决策边界的贡献为0？…该感谢上帝，还是感谢广大和伟大的科研工作者啊！让我等凡夫俗子可以瞥见如此精妙和无与伦比的数学之美！

       到这里，和支持向量机相关的东西就介绍完了。总结一下：支持向量机的基本思想可以概括为，首先通过非线性变换将输入空间变换到一个高维的空间，然后在这个新的空间求最优分类面即最大间隔分类面，而这种非线性变换是通过定义适当的内积核函数来实现的。SVM实际上是根据统计学习理论依照结构风险最小化的原则提出的，要求实现两个目的：1）两类问题能够分开（经验风险最小）2）margin最大化（风险上界最小）既是在保证风险最小的子集中选择经验风险最小的函数。

 

六、多类分类之SVM

       SVM是一种典型的两类分类器，即它只回答属于正类还是负类的问题。而现实中要解决的问题，往往是多类的问题。那如何由两类分类器得到多类分类器呢？

6.1、“一对多”的方法

       One-Against-All这个方法还是比较容易想到的。就是每次仍然解一个两类分类的问题。比如我们5个类别，第一次就把类别1的样本定为正样本，其余2，3，4，5的样本合起来定为负样本，这样得到一个两类分类器，它能够指出一个样本是还是不是第1类的；第二次我们把类别2 的样本定为正样本，把1，3，4，5的样本合起来定为负样本，得到一个分类器，如此下去，我们可以得到5个这样的两类分类器（总是和类别的数目一致）。到了有样本需要分类的时候，我们就拿着这个样本挨个分类器的问：是属于你的么？是属于你的么？哪个分类器点头说是了，文章的类别就确定了。这种方法的好处是每个优化问题的规模比较小，而且分类的时候速度很快（只需要调用5个分类器就知道了结果）。但有时也会出现两种很尴尬的情况，例如拿这个样本问了一圈，每一个分类器都说它是属于它那一类的，或者每一个分类器都说它不是它那一类的，前者叫分类重叠现象，后者叫不可分类现象。分类重叠倒还好办，随便选一个结果都不至于太离谱，或者看看这篇文章到各个超平面的距离，哪个远就判给哪个。不可分类现象就着实难办了，只能把它分给第6个类别了……更要命的是，本来各个类别的样本数目是差不多的，但“其余”的那一类样本数总是要数倍于正类（因为它是除正类以外其他类别的样本之和嘛），这就人为的造成了上一节所说的“数据集偏斜”问题。

       如下图左。红色分类面将红色与其他两种颜色分开，绿色分类面将绿色与其他两种颜色分开，蓝色分类面将蓝色与其他两种颜色分开。

       在这里的对某个点的分类实际上是通过衡量这个点到三个决策边界的距离，因为到分类面的距离越大，分类越可信嘛。当然了，这个距离是有符号的，如下所示：

       例如下图左，将星星这个点划分给绿色这一类。右图将星星这个点划分给褐色这一类。

6.2、“一对一”的方法

       One-Against-One方法是每次选一个类的样本作正类样本，而负类样本则变成只选一个类（称为“一对一单挑”的方法，哦，不对，没有单挑，就是“一对一”的方法，呵呵），这就避免了偏斜。因此过程就是算出这样一些分类器，第一个只回答“是第1类还是第2类”，第二个只回答“是第1类还是第3类”，第三个只回答“是第1类还是第4类”，如此下去，你也可以马上得出，这样的分类器应该有5 X 4/2=10个（通式是，如果有k个类别，则总的两类分类器数目为k(k-1)/2）。虽然分类器的数目多了，但是在训练阶段（也就是算出这些分类器的分类平面时）所用的总时间却比“一类对其余”方法少很多，在真正用来分类的时候，把一个样本扔给所有分类器，第一个分类器会投票说它是“1”或者“2”，第二个会说它是“1”或者“3”，让每一个都投上自己的一票，最后统计票数，如果类别“1”得票最多，就判这篇文章属于第1类。这种方法显然也会有分类重叠的现象，但不会有不可分类现象，因为总不可能所有类别的票数都是0。如下图右，中间紫色的块，每类的得票数都是1，那就不知道归类给那个类好了，只能随便扔给某个类了（或者衡量这个点到三个决策边界的距离，因为到分类面的距离越大，分类越可信嘛），扔掉了就是你命好，扔错了就不lucky了。

 

七、KKT条件分析

       对KKT条件，请大家参考文献[13][14]。假设我们优化得到的最优解是：α_i*,β_i*, ξ_i*, w*和b*。我们的最优解需要满足KKT条件：

      同时β_i*和ξ_i*都需要大于等于0，而α_i*需要在0和C之间。那可以分三种情况讨论：

       总的来说就是， KKT条件就变成了：

       第一个式子表明如果α_i=0，那么该样本落在两条间隔线外。第二个式子表明如果α_i=C，那么该样本有可能落在两条间隔线内部，也有可能落在两条间隔线上面，主要看对应的松弛变量的取值是等于0还是大于0，第三个式子表明如果0<α_i<C，那么该样本一定落在分隔线上（这点很重要，b就是拿这些落在分隔线上的点来求的，因为在分割线上w^Tx+b=1或者-1嘛，才是等式，在其他地方，都是不等式，求解不了b）。具体形象化的表示如下：

       通过KKT条件可知，α_i不等于0的都是支持向量，它有可能落在分隔线上，也有可能落在两条分隔线内部。KKT条件是非常重要的，在SMO也就是SVM的其中一个实现算法中，我们可以看到它的重要应用。

八、SVM的实现之SMO算法

      终于到SVM的实现部分了。那么神奇和有效的东西还得回归到实现才可以展示其强大的功力。SVM有效而且存在很高效的训练算法，这也是工业界非常青睐SVM的原因。

      前面讲到，SVM的学习问题可以转化为下面的对偶问题：

       需要满足的KKT条件：

       也就是说找到一组α_i可以满足上面的这些条件的就是该目标的一个最优解。所以我们的优化目标是找到一组最优的α_i*。一旦求出这些α_i*，就很容易计算出权重向量w*和b，并得到分隔超平面了。

       这是个凸二次规划问题，它具有全局最优解，一般可以通过现有的工具来优化。但当训练样本非常多的时候，这些优化算法往往非常耗时低效，以致无法使用。从SVM提出到现在，也出现了很多优化训练的方法。其中，非常出名的一个是1982年由Microsoft Research的John C. Platt在论文《Sequential Minimal Optimization: A Fast Algorithm for TrainingSupport Vector Machines》中提出的Sequential Minimal Optimization序列最小化优化算法，简称SMO算法。SMO算法的思想很简单，它将大优化的问题分解成多个小优化的问题。这些小问题往往比较容易求解，并且对他们进行顺序求解的结果与将他们作为整体来求解的结果完全一致。在结果完全一致的同时，SMO的求解时间短很多。在深入SMO算法之前，我们先来了解下坐标下降这个算法，SMO其实基于这种简单的思想的。

 

8.1、坐标下降（上升）法

      假设要求解下面的优化问题：

      在这里，我们需要求解m个变量α_i，一般来说是通过梯度下降（这里是求最大值，所以应该叫上升）等算法每一次迭代对所有m个变量α_i也就是α向量进行一次性优化。通过误差每次迭代调整α向量中每个元素的值。而坐标上升法（坐标上升与坐标下降可以看做是一对，坐标上升是用来求解max最优化问题，坐标下降用于求min最优化问题）的思想是每次迭代只调整一个变量α_i的值，其他变量的值在这次迭代中固定不变。

       最里面语句的意思是固定除α_i之外的所有α_j(i不等于j)，这时W可看作只是关于α_i的函数，那么直接对α_i求导优化即可。这里我们进行最大化求导的顺序i是从1到m，可以通过更改优化顺序来使W能够更快地增加并收敛。如果W在内循环中能够很快地达到最优，那么坐标上升法会是一个很高效的求极值方法。

      用个二维的例子来说明下坐标下降法：我们需要寻找f(x,y)=x²+xy+y²的最小值处的(x*, y*)，也就是下图的F*点的地方。

       假设我们初始的点是A（图是函数投影到xoy平面的等高线图，颜色越深值越小），我们需要达到F*的地方。那最快的方法就是图中黄色线的路径，一次性就到达了，其实这个是牛顿优化法，但如果是高维的话，这个方法就不太高效了（因为需要求解矩阵的逆，这个不在这里讨论）。我们也可以按照红色所指示的路径来走。从A开始，先固定x，沿着y轴往让f(x, y)值减小的方向走到B点，然后固定y，沿着x轴往让f(x, y)值减小的方向走到C点，不断循环，直到到达F*。反正每次只要我们都往让f(x, y)值小的地方走就行了，这样脚踏实地，一步步走，每一步都使f(x, y)慢慢变小，总有一天，皇天不负有心人的。到达F*也是时间问题。到这里你可能会说，这红色线比黄色线贫富差距也太严重了吧。因为这里是二维的简单的情况嘛。如果是高维的情况，而且目标函数很复杂的话，再加上样本集很多，那么在梯度下降中，目标函数对所有α_i求梯度或者在牛顿法中对矩阵求逆，都是很耗时的。这时候，如果W只对单个α_i优化很快的时候，坐标下降法可能会更加高效。

 

8.2、SMO算法

       SMO算法的思想和坐标下降法的思想差不多。唯一不同的是，SMO是一次迭代优化两个α而不是一个。为什么要优化两个呢？

       我们回到这个优化问题。我们可以看到这个优化问题存在着一个约束，也就是

       假设我们首先固定除α₁以外的所有参数，然后在α₁上求极值。但需要注意的是，因为如果固定α₁以外的所有参数，由上面这个约束条件可以知道，α₁将不再是变量（可以由其他值推出），因为问题中规定了：

      因此，我们需要一次选取两个参数做优化，比如α_i和α_j，此时α_i可以由α_j和其他参数表示出来。这样回代入W中，W就只是关于α_j的函数了，这时候就可以只对α_j进行优化了。在这里就是对α_j进行求导，令导数为0就可以解出这个时候最优的α_j了。然后也可以得到α_i。这就是一次的迭代过程，一次迭代只调整两个拉格朗日乘子α_i和α_j。SMO之所以高效就是因为在固定其他参数后，对一个参数优化过程很高效（对一个参数的优化可以通过解析求解，而不是迭代。虽然对一个参数的一次最小优化不可能保证其结果就是所优化的拉格朗日乘子的最终结果，但会使目标函数向极小值迈进一步，这样对所有的乘子做最小优化，直到所有满足KKT条件时，目标函数达到最小）。

       总结下来是：

重复下面过程直到收敛{

（1）选择两个拉格朗日乘子α_i和α_j；

（2）固定其他拉格朗日乘子α_k(k不等于i和j)，只对α_i和α_j优化w(α);

（3）根据优化后的α_i和α_j，更新截距b的值；

}

        那训练里面这两三步骤到底是怎么实现的，需要考虑什么呢？下面我们来具体分析下：

（1）选择α_i和α_j：

        我们现在是每次迭代都优化目标函数的两个拉格朗日乘子α_i和α_j，然后其他的拉格朗日乘子保持固定。如果有N个训练样本，我们就有N个拉格朗日乘子需要优化，但每次我们只挑两个进行优化，我们就有N(N-1)种选择。那到底我们要选择哪对α_i和α_j呢？选择哪对才好呢？想想我们的目标是什么？我们希望把所有违法KKT条件的样本都纠正回来，因为如果所有样本都满足KKT条件的话，我们的优化就完成了。那就很直观了，哪个害群之马最严重，我们得先对他进行思想教育，让他尽早回归正途。OK，那我们选择的第一个变量α_i就选违法KKT条件最严重的那一个。那第二个变量α_j怎么选呢？

       我们是希望快点找到最优的N个拉格朗日乘子，使得代价函数最大，换句话说，要最快的找到代价函数最大值的地方对应的N个拉格朗日乘子。这样我们的训练时间才会短。就像你从广州去北京，有飞机和绿皮车给你选，你选啥？（就算你不考虑速度，也得考虑下空姐的感受嘛，别辜负了她们渴望看到你的期盼，哈哈）。有点离题了，anyway，每次迭代中，哪对α_i和α_j可以让我更快的达到代价函数值最大的地方，我们就选他们。或者说，走完这一步，选这对α_i和α_j代价函数值增加的值最多，比选择其他所有α_i和α_j的结合中都多。这样我们才可以更快的接近代价函数的最大值，也就是达到优化的目标了。再例如，下图，我们要从A点走到B点，按蓝色的路线走c2方向的时候，一跨一大步，按红色的路线走c1方向的时候，只能是人类的一小步。所以，蓝色路线走两步就迈进了成功之门，而红色的路线，人生曲折，好像成功遥遥无期一样，故曰，选择比努力更重要！

       真啰嗦！说了半天，其实就一句话：为什么每次迭代都要选择最好的α_i和α_j，就是为了更快的收敛！那实践中每次迭代到底要怎样选α_i和α_j呢？这有个很好听的名字叫启发式选择，主要思想是先选择最有可能需要优化（也就是违反KKT条件最严重）的α_i，再针对这样的α_i选择最有可能取得较大修正步长的α_j。具体是以下两个过程：

1）第一个变量α_i的选择：

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

       该检验是在ε范围内进行的。在检验过程中，外层循环首先遍历所有满足条件0<α_j<C的样本点，即在间隔边界上的支持向量点，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的α_i。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的α_i。

       优先选择遍历非边界数据样本，因为非边界数据样本更有可能需要调整，边界数据样本常常不能得到进一步调整而留在边界上。由于大部分数据样本都很明显不可能是支持向量，因此对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本得到调整时，遍历所有数据样本，以检验是否整个集合都满足KKT条件。如果整个集合的检验中又有数据样本被进一步进化，则有必要再遍历非边界数据样本。这样，不停地在遍历所有数据样本和遍历非边界数据样本之间切换，直到整个样本集合都满足KKT条件为止。以上用KKT条件对数据样本所做的检验都以达到一定精度ε就可以停止为条件。如果要求十分精确的输出算法，则往往不能很快收敛。

       对整个数据集的遍历扫描相当容易，而实现对非边界α_i的扫描时，首先需要将所有非边界样本的α_i值（也就是满足0<α_i<C）保存到新的一个列表中，然后再对其进行遍历。同时，该步骤跳过那些已知的不会改变的α_i值。

2）第二个变量α_j的选择：

       在选择第一个α_i后，算法会通过一个内循环来选择第二个α_j值。因为第二个乘子的迭代步长大致正比于|E_i-E_j|，所以我们需要选择能够最大化|E_i-E_j|的第二个乘子（选择最大化迭代步长的第二个乘子）。在这里，为了节省计算时间，我们建立一个全局的缓存用于保存所有样本的误差值，而不用每次选择的时候就重新计算。我们从中选择使得步长最大或者|E_i-E_j|最大的α_j。

（2）优化α_i和α_j：

       选择这两个拉格朗日乘子后，我们需要先计算这些参数的约束值。然后再求解这个约束最大化问题。

       首先，我们需要给α_j找到边界L<=α_j<=H，以保证α_j满足0<=α_j<=C的约束。这意味着α_j必须落入这个盒子中。由于只有两个变量(α_i, α_j)，约束可以用二维空间中的图形来表示，如下图：

       不等式约束使得(α_i,α_j)在盒子[0, C]x[0, C]内，等式约束使得(α_i, α_j)在平行于盒子[0, C]x[0, C]的对角线的直线上。因此要求的是目标函数在一条平行于对角线的线段上的最优值。这使得两个变量的最优化问题成为实质的单变量的最优化问题。由图可以得到，α_j的上下界可以通过下面的方法得到：

       我们优化的时候，α_j必须要满足上面这个约束。也就是说上面是α_j的可行域。然后我们开始寻找α_j，使得目标函数最大化。通过推导得到α_j的更新公式如下：

       这里E_k可以看做对第k个样本，SVM的输出与期待输出，也就是样本标签的误差。

       而η实际上是度量两个样本i和j的相似性的。在计算η的时候，我们需要使用核函数，那么就可以用核函数来取代上面的内积。

       得到新的α_j后，我们需要保证它处于边界内。换句话说，如果这个优化后的值跑出了边界L和H，我们就需要简单的裁剪，将α_j收回这个范围：

       最后，得到优化的α_j后，我们需要用它来计算α_i：

       到这里，α_i和α_j的优化就完成了。

（3）计算阈值b：

       优化α_i和α_j后，我们就可以更新阈值b，使得对两个样本i和j都满足KKT条件。如果优化后α_i不在边界上（也就是满足0<α_i<C，这时候根据KKT条件，可以得到y_ig_i(x_i)=1，这样我们才可以计算b），那下面的阈值b1是有效的，因为当输入x_i时它迫使SVM输出y_i。

       同样，如果0<α_j<C，那么下面的b2也是有效的：

      如果0<α_i<C和0<α_j<C都满足，那么b1和b2都有效，而且他们是相等的。如果他们两个都处于边界上（也就是α_i=0或者α_i=C，同时α_j=0或者α_j=C），那么在b1和b2之间的阈值都满足KKT条件，一般我们取他们的平均值b=(b1+b2)/2。所以，总的来说对b的更新如下：

       每做完一次最小优化，必须更新每个数据样本的误差，以便用修正过的分类面对其他数据样本再做检验，在选择第二个配对优化数据样本时用来估计步长。

（4）凸优化问题终止条件：

       SMO算法的基本思路是：如果说有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件（证明请参考文献）。所以我们可以监视原问题的KKT条件，所以所有的样本都满足KKT条件，那么就表示迭代结束了。但是由于KKT条件本身是比较苛刻的，所以也需要设定一个容忍值，即所有样本在容忍值范围内满足KKT条件则认为训练可以结束；当然了，对于对偶问题的凸优化还有其他终止条件，可以参考文献。

 

8.3、SMO算法的Python实现

8.3.1、Python的准备工作

      我使用的Python是2.7.5版本的。附加的库有Numpy和Matplotlib。而Matplotlib又依赖dateutil和pyparsing两个库，所以我们需要安装以上三个库。前面两个库还好安装，直接在官网下对应版本就行。但我找后两个库的时候，就没那么容易了。后来发现，其实对Python的库的下载和安装可以借助pip工具的。这个是安装和管理Python包的工具。感觉它有点像ubuntu的apt-get，需要安装什么库，直接下载和安装一条龙服务。

       首先，我们需要到pip的官网：https://pypi.python.org/pypi/pip下载对应我们python版本的pip，例如我的是pip-1.4.1.tar.gz。但安装pip需要另一个工具，也就是setuptools，我们到https://pypi.python.org/pypi/setuptools/#windows下载ez_setup.py这个文件回来。然后在CMD命令行中执行：（注意他们的路径）

#python ez_setup.py

这时候，就会自动下载.egg等等文件然后安装完成。

      然后我们解压pip-1.4.1.tar.gz。进入到该目录中，执行：

#python setup.py install

这时候就会自动安装pip到你python目录下的Scripts文件夹中。我的是C:\Python27\Scripts。

       在里面我们可以看到pip.exe，然后我们进入到该文件夹中：

#cd C:\Python27\Scripts

#pip install dateutil

#pip install pyparsing

这样就可以把这些额外的库给下载回来了。非常高端大气上档次！

8.3.2、SMO算法的Python实现

       在代码中已经有了比较详细的注释了。不知道有没有错误的地方，如果有，还望大家指正（每次的运行结果都有可能不同，另外，感觉有些结果似乎不太正确，但我还没发现哪里出错了，如果大家找到有错误的地方，还望大家指点下，衷心感谢）。里面我写了个可视化结果的函数，但只能在二维的数据上面使用。直接贴代码：

SVM.py

[python] view plaincopy

1. #################################################  
2. # SVM: support vector machine  
3. # Author : zouxy  
4. # Date   : 2013-12-12  
5. # HomePage : http://blog.csdn.net/zouxy09  
6. # Email  : [e-mail]  
7. #################################################  
8.   
9. from numpy import *  
10. import time  
11. import matplotlib.pyplot as plt   
12.   
13.   
14. # calulate kernel value  
15. def calcKernelValue(matrix_x, sample_x, kernelOption):  
16.     kernelType = kernelOption[0]  
17.     numSamples = matrix_x.shape[0]  
18.     kernelValue = mat(zeros((numSamples, 1)))  
19.       
20.     if kernelType == 'linear':  
21.         kernelValue = matrix_x * sample_x.T  
22.     elif kernelType == 'rbf':  
23.         sigma = kernelOption[1]  
24.         if sigma == 0:  
25.             sigma = 1.0  
26.         for i in xrange(numSamples):  
27.             diff = matrix_x[i, :] - sample_x  
28.             kernelValue[i] = exp(diff * diff.T / (-2.0 * sigma**2))  
29.     else:  
30.         raise NameError('Not support kernel type! You can use linear or rbf!')  
31.     return kernelValue  
32.   
33.   
34. # calculate kernel matrix given train set and kernel type  
35. def calcKernelMatrix(train_x, kernelOption):  
36.     numSamples = train_x.shape[0]  
37.     kernelMatrix = mat(zeros((numSamples, numSamples)))  
38.     for i in xrange(numSamples):  
39.         kernelMatrix[:, i] = calcKernelValue(train_x, train_x[i, :], kernelOption)  
40.     return kernelMatrix  
41.   
42.   
43. # define a struct just for storing variables and data  
44. class SVMStruct:  
45.     def __init__(self, dataSet, labels, C, toler, kernelOption):  
46.         self.train_x = dataSet # each row stands for a sample  
47.         self.train_y = labels  # corresponding label  
48.         self.C = C             # slack variable  
49.         self.toler = toler     # termination condition for iteration  
50.         self.numSamples = dataSet.shape[0] # number of samples  
51.         self.alphas = mat(zeros((self.numSamples, 1))) # Lagrange factors for all samples  
52.         self.b = 0  
53.         self.errorCache = mat(zeros((self.numSamples, 2)))  
54.         self.kernelOpt = kernelOption  
55.         self.kernelMat = calcKernelMatrix(self.train_x, self.kernelOpt)  
56.   
57.           
58. # calculate the error for alpha k  
59. def calcError(svm, alpha_k):  
60.     output_k = float(multiply(svm.alphas, svm.train_y).T * svm.kernelMat[:, alpha_k] + svm.b)  
61.     error_k = output_k - float(svm.train_y[alpha_k])  
62.     return error_k  
63.   
64.   
65. # update the error cache for alpha k after optimize alpha k  
66. def updateError(svm, alpha_k):  
67.     error = calcError(svm, alpha_k)  
68.     svm.errorCache[alpha_k] = [1, error]  
69.   
70.   
71. # select alpha j which has the biggest step  
72. def selectAlpha_j(svm, alpha_i, error_i):  
73.     svm.errorCache[alpha_i] = [1, error_i] # mark as valid(has been optimized)  
74.     candidateAlphaList = nonzero(svm.errorCache[:, 0].A)[0] # mat.A return array  
75.     maxStep = 0; alpha_j = 0; error_j = 0  
76.   
77.     # find the alpha with max iterative step  
78.     if len(candidateAlphaList) > 1:  
79.         for alpha_k in candidateAlphaList:  
80.             if alpha_k == alpha_i:   
81.                 continue  
82.             error_k = calcError(svm, alpha_k)  
83.             if abs(error_k - error_i) > maxStep:  
84.                 maxStep = abs(error_k - error_i)  
85.                 alpha_j = alpha_k  
86.                 error_j = error_k  
87.     # if came in this loop first time, we select alpha j randomly  
88.     else:             
89.         alpha_j = alpha_i  
90.         while alpha_j == alpha_i:  
91.             alpha_j = int(random.uniform(0, svm.numSamples))  
92.         error_j = calcError(svm, alpha_j)  
93.       
94.     return alpha_j, error_j  
95.   
96.   
97. # the inner loop for optimizing alpha i and alpha j  
98. def innerLoop(svm, alpha_i):  
99.     error_i = calcError(svm, alpha_i)  
100.   
101.     ### check and pick up the alpha who violates the KKT condition  
102.     ## satisfy KKT condition  
103.     # 1) yi*f(i) >= 1 and alpha == 0 (outside the boundary)  
104.     # 2) yi*f(i) == 1 and 0<alpha< C (on the boundary)  
105.     # 3) yi*f(i) <= 1 and alpha == C (between the boundary)  
106.     ## violate KKT condition  
107.     # because y[i]*E_i = y[i]*f(i) - y[i]^2 = y[i]*f(i) - 1, so  
108.     # 1) if y[i]*E_i < 0, so yi*f(i) < 1, if alpha < C, violate!(alpha = C will be correct)   
109.     # 2) if y[i]*E_i > 0, so yi*f(i) > 1, if alpha > 0, violate!(alpha = 0 will be correct)  
110.     # 3) if y[i]*E_i = 0, so yi*f(i) = 1, it is on the boundary, needless optimized  
111.     if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\  
112.         (svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):  
113.   
114.         # step 1: select alpha j  
115.         alpha_j, error_j = selectAlpha_j(svm, alpha_i, error_i)  
116.         alpha_i_old = svm.alphas[alpha_i].copy()  
117.         alpha_j_old = svm.alphas[alpha_j].copy()  
118.   
119.         # step 2: calculate the boundary L and H for alpha j  
120.         if svm.train_y[alpha_i] != svm.train_y[alpha_j]:  
121.             L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])  
122.             H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])  
123.         else:  
124.             L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)  
125.             H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])  
126.         if L == H:  
127.             return 0  
128.   
129.         # step 3: calculate eta (the similarity of sample i and j)  
130.         eta = 2.0 * svm.kernelMat[alpha_i, alpha_j] - svm.kernelMat[alpha_i, alpha_i] \  
131.                   - svm.kernelMat[alpha_j, alpha_j]  
132.         if eta >= 0:  
133.             return 0  
134.   
135.         # step 4: update alpha j  
136.         svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta  
137.   
138.         # step 5: clip alpha j  
139.         if svm.alphas[alpha_j] > H:  
140.             svm.alphas[alpha_j] = H  
141.         if svm.alphas[alpha_j] < L:  
142.             svm.alphas[alpha_j] = L  
143.   
144.         # step 6: if alpha j not moving enough, just return       
145.         if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:  
146.             updateError(svm, alpha_j)  
147.             return 0  
148.   
149.         # step 7: update alpha i after optimizing aipha j  
150.         svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \  
151.                                 * (alpha_j_old - svm.alphas[alpha_j])  
152.   
153.         # step 8: update threshold b  
154.         b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \  
155.                                                     * svm.kernelMat[alpha_i, alpha_i] \  
156.                              - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \  
157.                                                     * svm.kernelMat[alpha_i, alpha_j]  
158.         b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \  
159.                                                     * svm.kernelMat[alpha_i, alpha_j] \  
160.                              - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \  
161.                                                     * svm.kernelMat[alpha_j, alpha_j]  
162.         if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):  
163.             svm.b = b1  
164.         elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):  
165.             svm.b = b2  
166.         else:  
167.             svm.b = (b1 + b2) / 2.0  
168.   
169.         # step 9: update error cache for alpha i, j after optimize alpha i, j and b  
170.         updateError(svm, alpha_j)  
171.         updateError(svm, alpha_i)  
172.   
173.         return 1  
174.     else:  
175.         return 0  
176.   
177.   
178. # the main training procedure  
179. def trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('rbf', 1.0)):  
180.     # calculate training time  
181.     startTime = time.time()  
182.   
183.     # init data struct for svm  
184.     svm = SVMStruct(mat(train_x), mat(train_y), C, toler, kernelOption)  
185.       
186.     # start training  
187.     entireSet = True  
188.     alphaPairsChanged = 0  
189.     iterCount = 0  
190.     # Iteration termination condition:  
191.     #   Condition 1: reach max iteration  
192.     #   Condition 2: no alpha changed after going through all samples,  
193.     #                in other words, all alpha (samples) fit KKT condition  
194.     while (iterCount < maxIter) and ((alphaPairsChanged > 0) or entireSet):  
195.         alphaPairsChanged = 0  
196.   
197.         # update alphas over all training examples  
198.         if entireSet:  
199.             for i in xrange(svm.numSamples):  
200.                 alphaPairsChanged += innerLoop(svm, i)  
201.             print '---iter:%d entire set, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)  
202.             iterCount += 1  
203.         # update alphas over examples where alpha is not 0 & not C (not on boundary)  
204.         else:  
205.             nonBoundAlphasList = nonzero((svm.alphas.A > 0) * (svm.alphas.A < svm.C))[0]  
206.             for i in nonBoundAlphasList:  
207.                 alphaPairsChanged += innerLoop(svm, i)  
208.             print '---iter:%d non boundary, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)  
209.             iterCount += 1  
210.   
211.         # alternate loop over all examples and non-boundary examples  
212.         if entireSet:  
213.             entireSet = False  
214.         elif alphaPairsChanged == 0:  
215.             entireSet = True  
216.   
217.     print 'Congratulations, training complete! Took %fs!' % (time.time() - startTime)  
218.     return svm  
219.   
220.   
221. # testing your trained svm model given test set  
222. def testSVM(svm, test_x, test_y):  
223.     test_x = mat(test_x)  
224.     test_y = mat(test_y)  
225.     numTestSamples = test_x.shape[0]  
226.     supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]  
227.     supportVectors      = svm.train_x[supportVectorsIndex]  
228.     supportVectorLabels = svm.train_y[supportVectorsIndex]  
229.     supportVectorAlphas = svm.alphas[supportVectorsIndex]  
230.     matchCount = 0  
231.     for i in xrange(numTestSamples):  
232.         kernelValue = calcKernelValue(supportVectors, test_x[i, :], svm.kernelOpt)  
233.         predict = kernelValue.T * multiply(supportVectorLabels, supportVectorAlphas) + svm.b  
234.         if sign(predict) == sign(test_y[i]):  
235.             matchCount += 1  
236.     accuracy = float(matchCount) / numTestSamples  
237.     return accuracy  
238.   
239.   
240. # show your trained svm model only available with 2-D data  
241. def showSVM(svm):  
242.     if svm.train_x.shape[1] != 2:  
243.         print "Sorry! I can not draw because the dimension of your data is not 2!"  
244.         return 1  
245.   
246.     # draw all samples  
247.     for i in xrange(svm.numSamples):  
248.         if svm.train_y[i] == -1:  
249.             plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'or')  
250.         elif svm.train_y[i] == 1:  
251.             plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'ob')  
252.   
253.     # mark support vectors  
254.     supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]  
255.     for i in supportVectorsIndex:  
256.         plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'oy')  
257.       
258.     # draw the classify line  
259.     w = zeros((2, 1))  
260.     for i in supportVectorsIndex:  
261.         w += multiply(svm.alphas[i] * svm.train_y[i], svm.train_x[i, :].T)   
262.     min_x = min(svm.train_x[:, 0])[0, 0]  
263.     max_x = max(svm.train_x[:, 0])[0, 0]  
264.     y_min_x = float(-svm.b - w[0] * min_x) / w[1]  
265.     y_max_x = float(-svm.b - w[0] * max_x) / w[1]  
266.     plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')  
267.     plt.show()  

       测试的数据来自这里。有100个样本，每个样本两维，最后是对应的标签，例如：

3.542485 1.977398          -1

3.018896 2.556416          -1

7.551510 -1.580030         1

2.114999 -0.004466         -1

……

       测试代码中首先加载这个数据库，然后用前面80个样本来训练，再用剩下的20个样本的测试，并显示训练后的模型和分类结果。测试代码如下：

test_SVM.py

[python] view plaincopy

1. #################################################  
2. # SVM: support vector machine  
3. # Author : zouxy  
4. # Date   : 2013-12-12  
5. # HomePage : http://blog.csdn.net/zouxy09  
6. # Email  : [e-mail]  
7. #################################################  
8.   
9. from numpy import *  
10. import SVM  
11.   
12. ################## test svm #####################  
13. ## step 1: load data  
14. print "step 1: load data..."  
15. dataSet = []  
16. labels = []  
17. fileIn = open('E:/Python/Machine Learning in Action/testSet.txt')  
18. for line in fileIn.readlines():  
19.     lineArr = line.strip().split('\t')  
20.     dataSet.append([float(lineArr[0]), float(lineArr[1])])  
21.     labels.append(float(lineArr[2]))  
22.   
23. dataSet = mat(dataSet)  
24. labels = mat(labels).T  
25. train_x = dataSet[0:81, :]  
26. train_y = labels[0:81, :]  
27. test_x = dataSet[80:101, :]  
28. test_y = labels[80:101, :]  
29.   
30. ## step 2: training...  
31. print "step 2: training..."  
32. C = 0.6  
33. toler = 0.001  
34. maxIter = 50  
35. svmClassifier = SVM.trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('linear', 0))  
36.   
37. ## step 3: testing  
38. print "step 3: testing..."  
39. accuracy = SVM.testSVM(svmClassifier, test_x, test_y)  
40.   
41. ## step 4: show the result  
42. print "step 4: show the result..."    
43. print 'The classify accuracy is: %.3f%%' % (accuracy * 100)  
44. SVM.showSVM(svmClassifier)  

运行结果如下：

[python] view plaincopy

1. step 1: load data...  
2. step 2: training...  
3. ---iter:0 entire set, alpha pairs changed:8  
4. ---iter:1 non boundary, alpha pairs changed:7  
5. ---iter:2 non boundary, alpha pairs changed:1  
6. ---iter:3 non boundary, alpha pairs changed:0  
7. ---iter:4 entire set, alpha pairs changed:0  
8. Congratulations, training complete! Took 0.058000s!  
9. step 3: testing...  
10. step 4: show the result...  
11. The classify accuracy is: 100.000%  

训练好的模型图：

 

 

九、参考文献与推荐阅读

[1] JerryLead的博客，作者根据斯坦福的讲义给出了流畅和通俗的推导：SVM系列。

[2]嘉士伯的SVM入门系列，讲得很好。

[3] pluskid的支持向量机系列，非常好。其中关于dual问题推导非常赞。

[4] Leo Zhang的SVM学习系列，博客中还包含了很多其他的机器学习算法。

[5] v_july_v的支持向量机通俗导论（理解SVM的三层境界）。结构之法算法之道blog。

[6] 李航的《统计学习方法》，清华大学出版社

[7] SVM学习——Sequential Minimal Optimization

[8] SVM算法实现（一）

[9] Sequential Minimal Optimization: A FastAlgorithm for Training Support Vector Machines

[10] SVM --从“原理”到实现

[11] 支持向量机入门系列

[12]SVM的各个版本及其多种语言实现代码合集

[13] Karush-Kuhn-Tucker(KKT) conditions

[14] 深入理解拉格朗日乘子法（Lagrange Multiplier) 和KKT条件