机器学习之支持向量机—SVM原理代码实现
支持向量机—SVM原理代码实现
本文系作者原创，转载请注明出处:https://www.cnblogs.com/further-further-further/p/9596898.html
1. 解决什么问题？
最基本的应用是数据分类，特别是对于非线性不可分数据集。支持向量机不仅能对非线性可分数据集进行分类，对于非线性不可分数据集的也可以分类
（我认为这才是支持向量机的真正魅力所在，因为现实场景中，样本数据往往是非线性不可分的）。
现实场景一 ：样本数据大部分是线性可分的，但是只是在样本中含有少量噪声或特异点，去掉这些噪声或特异点后线性可分 => 用支持向量机的软间隔方法进行分类；
现实场景二：样本数据完全线性不可分 => 引入核函数，将低维不可分的非线性数据集转化为高维可分的数据集，用支持向量机的软间隔方法进行分类；
一定要明白一点：对分类器准确性有影响只是样本中的支持向量，因此，其他样本在计算分类器的权重矩阵时可以直接过滤掉，大大节省运行时间。
目标：获取离边界最近样本点到超平面最远距离
2. SVM的思想演化？
博客 https://www.cnblogs.com/zhizhan/p/4430253.html 关于SVM的演化说得很透彻，也很形象，这里借用一下。
2.1 硬间隔支持向量机
SVM中最关键的思想之一就是引入和定义了“间隔”这个概念。这个概念本身很简单，以二维空间为例，就是点到分类直线之间的距离。假设直线为y=wx+b，那么只要使所有正分类点到该直线的距离与所有负分类点到该直线的距离的总和达到最大，这条直线就是最优分类直线。这样，原问题就转化为一个约束优化问题，可以直接求解。这叫做硬间隔最大化，得到的SVM模型称作硬间隔支持向量机。
2.2 软间隔支持向量机
但是新问题出现了，在实际应用中，我们得到的数据并不总是完美的线性可分的，其中可能会有个别噪声点，他们错误的被分类到了其他类中。如果将这些特异的噪点去除后，可以很容易的线性可分。但是，我们对于数据集中哪些是噪声点却是不知道的，如果以之前的方法进行求解，会无法进行线性分开。是不是就没办法了呢？假设在y=x+1直线上下分为两类，若两类中各有对方的几个噪点，在人的眼中，仍然是可以将两类分开的。这是因为在人脑中是可以容忍一定的误差的，仍然使用y=x+1直线分类，可以在最小误差的情况下进行最优的分类。同样的道理，我们在SVM中引入误差的概念，将其称作“松弛变量”。通过加入松弛变量，在原距离函数中需要加入新的松弛变量带来的误差，这样，最终的优化目标函数变成了两个部分组成：距离函数和松弛变量误差。这两个部分的重要程度并不是相等的，而是需要依据具体问题而定的，因此，我们加入权重参数C，将其与目标函数中的松弛变量误差相乘，这样，就可以通过调整C来对二者的系数进行调和。如果我们能够容忍噪声，那就把C调小，让他的权重降下来，从而变得不重要；反之，我们需要很严格的噪声小的模型，则将C调大一点，权重提升上去，变得更加重要。通过对参数C的调整，可以对模型进行控制。这叫做软间隔最大化，得到的SVM称作软间隔支持向量机。
2.3 非线性支持向量机
之前的硬间隔支持向量机和软间隔支持向量机都是解决线性可分数据集或近似线性可分数据集的问题的。但是如果噪点很多，甚至会造成数据变成了线性不可分的，那该怎么办？最常见的例子是在二维平面笛卡尔坐标系下，以原点(0,0)为圆心，以1为半径画圆，则圆内的点和圆外的点在二维空间中是肯定无法线性分开的。但是，学过初中几何就知道，对于圆圈内（含圆圈）的点：x^2+y^2≤1，圆圈外的则x^2+y^2＞1。我们假设第三个维度：z=x^2+y^2，那么在第三维空间中，可以通过z是否大于1来判断该点是否在圆内还是圆外。这样，在二维空间中线性不可分的数据在第三维空间很容易的线性可分了。这就是非线性支持向量机。
实际中，对某个实际问题函数来寻找一个合适的空间进行映射是非常困难的，幸运的是，在计算中发现，我们需要的只是两个向量在新的映射空间中的内积结果，而映射函数到底是怎么样的其实并不需要知道。这一点不太好理解，有人会问，既然不知道映射函数，那怎么能知道映射后在新空间中的内积结果呢？答案其实是可以的。这就需要引入了核函数的概念。核函数是这样的一种函数：仍然以二维空间为例，假设对于变量x和y，将其映射到新空间的映射函数为φ，则在新空间中，二者分别对应φ(x)和φ(y)，他们的内积则为<φ(x),φ(y)>。我们令函数Kernel(x,y)=<φ(x),φ(y)>=k(x,y)，可以看出，函数Kernel(x,y)是一个关于x和y的函数！而与φ无关！这是一个多么好的性质！我们再也不用管φ具体是什么映射关系了，只需要最后计算Kernel(x,y)就可以得到他们在高维空间中的内积，这样就可以直接带入之前的支持向量机中计算！
核函数不是很好找到，一般是由数学家反向推导出来或拼凑出来的。现在知道的线性核函数有多项式核函数、高斯核函数等。其中，高斯核函数对应的支持向量机是高斯径向基函数（RBF），是最常用的核函数。
RBF核函数可以将维度扩展到无穷维的空间，因此，理论上讲可以满足一切映射的需求。为什么会是无穷维呢？我以前都不太明白这一点。后来老师讲到，RBF对应的是泰勒级数展开，在泰勒级数中，一个函数可以分解为无穷多个项的加和，其中，每一个项可以看做是对应的一个维度，这样，原函数就可以看做是映射到了无穷维的空间中。这样，在实际应用中，RBF是相对最好的一个选择。当然，如果有研究的话，还可以选用其他核函数，可能会在某些问题上表现更好。但是，RBF是在对问题不了解的情况下，对最广泛问题效果都很不错的核函数。因此，使用范围也最广。
3. SVM原理
3.1 目标函数
求使几何间隔最大的分离超平面
以及相应的分类决策函数
也就是求出w，b
3.2 几何间隔定义
3.3 学习的对偶算法
为了求解线性可分支持向量机的最优化问题，将它作为原始最优化问题，应用拉格朗日对偶性，通过求解对偶问题得到原始问题最优解。
3.4 最大软间隔分类器（线性可分支持向量学习的最优化问题）
这就是线性可分支持向量机的对偶算法，这样做的优点，一是对偶问题往往更容易求解；二是自然引入核函数，进而推广到非线性分类问题。
的求解就转化为
的求解。上述公式就是对不同变量求偏导，具体推导过程见《统计学习方法》第7章 或者相关博客。
3.5 KKT条件
上图的总结有点小问题，
分类不能合并，正确的分类描述应该是：
在超平面上或误分的样本点是不能正确分类的。
4. SMO算法（sequential minimal optimization 序列最小最优化）
高效实现支持向量机的算法是SMO算法，其基本思路是：如果所有变量的解都满足此最优化的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该优化问题的充分必要条件。否则，
选择两个变量，固定其他变量，针对这两个变量构建一个二次规划问题，这个二次规划问题关于这两个变量的解应该更接近原始二次规划问题的解，因为这会使原始二次规划问题的目标函数值变得更小。
重要的是，这时子问题可以通过解析方法求解，这样可以大大提高整个算法的计算速度。
子问题有两个变量，一个是违反KKT条件最严重的那一个，另一个由约束条件自动确定。如此，SMO算法将原问题不断分解为子问题并对子问题求解，进而达到求解原问题的目的。
由
并且 推出
4.1 
范围
由
假定不固定，其他拉格朗日因子固定（也就是常量），得到
常量 ，
为变量，
为因变量，L，H的范围就是
的 范围
4.2 计算拉格朗日乘子 
，
第一个 公式如下，推导过程见《统计学习方法》第7章或者相关博客。
选择一个违背KKT的数据项，根据下面结论得到优化后 ，
s表示
，从而可以得到另一个优化后的 ，同样也需要进行约束。
4.3 计算b1，b2
在每次完成两个变量的优化后，都要重新计算阈值b。具体推导见《统计学习方法》
5. 代码实现
5.1 输入数据
这里使用两个原始数据文件 trainingData.txt，testData.txt。
trainingData.txt
1 -0.214824 0.662756 -1.000000 2 -0.061569 -0.091875 1.000000 3 0.406933 0.648055 -1.000000 4 0.223650 0.130142 1.000000 5 0.231317 0.766906 -1.000000 6 -0.748800 -0.531637 -1.000000 7 -0.557789 0.375797 -1.000000 8 0.207123 -0.019463 1.000000 9 0.286462 0.719470 -1.000000 10 0.195300 -0.179039 1.000000 11 -0.152696 -0.153030 1.000000 12 0.384471 0.653336 -1.000000 13 -0.117280 -0.153217 1.000000 14 -0.238076 0.000583 1.000000 15 -0.413576 0.145681 1.000000 16 0.490767 -0.680029 -1.000000 17 0.199894 -0.199381 1.000000 18 -0.356048 0.537960 -1.000000 19 -0.392868 -0.125261 1.000000 20 0.353588 -0.070617 1.000000 21 0.020984 0.925720 -1.000000 22 -0.475167 -0.346247 -1.000000 23 0.074952 0.042783 1.000000 24 0.394164 -0.058217 1.000000 25 0.663418 0.436525 -1.000000 26 0.402158 0.577744 -1.000000 27 -0.449349 -0.038074 1.000000 28 0.619080 -0.088188 -1.000000 29 0.268066 -0.071621 1.000000 30 -0.015165 0.359326 1.000000 31 0.539368 -0.374972 -1.000000 32 -0.319153 0.629673 -1.000000 33 0.694424 0.641180 -1.000000 34 0.079522 0.193198 1.000000 35 0.253289 -0.285861 1.000000 36 -0.035558 -0.010086 1.000000 37 -0.403483 0.474466 -1.000000 38 -0.034312 0.995685 -1.000000 39 -0.590657 0.438051 -1.000000 40 -0.098871 -0.023953 1.000000 41 -0.250001 0.141621 1.000000 42 -0.012998 0.525985 -1.000000 43 0.153738 0.491531 -1.000000 44 0.388215 -0.656567 -1.000000 45 0.049008 0.013499 1.000000 46 0.068286 0.392741 1.000000 47 0.747800 -0.066630 -1.000000 48 0.004621 -0.042932 1.000000 49 -0.701600 0.190983 -1.000000 50 0.055413 -0.024380 1.000000 51 0.035398 -0.333682 1.000000 52 0.211795 0.024689 1.000000 53 -0.045677 0.172907 1.000000 54 0.595222 0.209570 -1.000000 55 0.229465 0.250409 1.000000 56 -0.089293 0.068198 1.000000 57 0.384300 -0.176570 1.000000 58 0.834912 -0.110321 -1.000000 59 -0.307768 0.503038 -1.000000 60 -0.777063 -0.348066 -1.000000 61 0.017390 0.152441 1.000000 62 -0.293382 -0.139778 1.000000 63 -0.203272 0.286855 1.000000 64 0.957812 -0.152444 -1.000000 65 0.004609 -0.070617 1.000000 66 -0.755431 0.096711 -1.000000 67 -0.526487 0.547282 -1.000000 68 -0.246873 0.833713 -1.000000 69 0.185639 -0.066162 1.000000 70 0.851934 0.456603 -1.000000 71 -0.827912 0.117122 -1.000000 72 0.233512 -0.106274 1.000000 73 0.583671 -0.709033 -1.000000 74 -0.487023 0.625140 -1.000000 75 -0.448939 0.176725 1.000000 76 0.155907 -0.166371 1.000000 77 0.334204 0.381237 -1.000000 78 0.081536 -0.106212 1.000000 79 0.227222 0.527437 -1.000000 80 0.759290 0.330720 -1.000000 81 0.204177 -0.023516 1.000000 82 0.577939 0.403784 -1.000000 83 -0.568534 0.442948 -1.000000 84 -0.011520 0.021165 1.000000 85 0.875720 0.422476 -1.000000 86 0.297885 -0.632874 -1.000000 87 -0.015821 0.031226 1.000000 88 0.541359 -0.205969 -1.000000 89 -0.689946 -0.508674 -1.000000 90 -0.343049 0.841653 -1.000000 91 0.523902 -0.436156 -1.000000 92 0.249281 -0.711840 -1.000000 93 0.193449 0.574598 -1.000000 94 -0.257542 -0.753885 -1.000000 95 -0.021605 0.158080 1.000000 96 0.601559 -0.727041 -1.000000 97 -0.791603 0.095651 -1.000000 98 -0.908298 -0.053376 -1.000000 99 0.122020 0.850966 -1.000000 100 -0.725568 -0.292022 -1.000000
testData.txt
1 3.542485 1.977398 -1 2 3.018896 2.556416 -1 3 7.551510 -1.580030 1 4 2.114999 -0.004466 -1 5 8.127113 1.274372 1 6 7.108772 -0.986906 1 7 8.610639 2.046708 1 8 2.326297 0.265213 -1 9 3.634009 1.730537 -1 10 0.341367 -0.894998 -1 11 3.125951 0.293251 -1 12 2.123252 -0.783563 -1 13 0.887835 -2.797792 -1 14 7.139979 -2.329896 1 15 1.696414 -1.212496 -1 16 8.117032 0.623493 1 17 8.497162 -0.266649 1 18 4.658191 3.507396 -1 19 8.197181 1.545132 1 20 1.208047 0.213100 -1 21 1.928486 -0.321870 -1 22 2.175808 -0.014527 -1 23 7.886608 0.461755 1 24 3.223038 -0.552392 -1 25 3.628502 2.190585 -1 26 7.407860 -0.121961 1 27 7.286357 0.251077 1 28 2.301095 -0.533988 -1 29 -0.232542 -0.547690 -1 30 3.457096 -0.082216 -1 31 3.023938 -0.057392 -1 32 8.015003 0.885325 1 33 8.991748 0.923154 1 34 7.916831 -1.781735 1 35 7.616862 -0.217958 1 36 2.450939 0.744967 -1 37 7.270337 -2.507834 1 38 1.749721 -0.961902 -1 39 1.803111 -0.176349 -1 40 8.804461 3.044301 1 41 1.231257 -0.568573 -1 42 2.074915 1.410550 -1 43 -0.743036 -1.736103 -1 44 3.536555 3.964960 -1 45 8.410143 0.025606 1 46 7.382988 -0.478764 1 47 6.960661 -0.245353 1 48 8.234460 0.701868 1 49 8.168618 -0.903835 1 50 1.534187 -0.622492 -1 51 9.229518 2.066088 1 52 7.886242 0.191813 1 53 2.893743 -1.643468 -1 54 1.870457 -1.040420 -1 55 5.286862 -2.358286 1 56 6.080573 0.418886 1 57 2.544314 1.714165 -1 58 6.016004 -3.753712 1 59 0.926310 -0.564359 -1 60 0.870296 -0.109952 -1 61 2.369345 1.375695 -1 62 1.363782 -0.254082 -1 63 7.279460 -0.189572 1 64 1.896005 0.515080 -1 65 8.102154 -0.603875 1 66 2.529893 0.662657 -1 67 1.963874 -0.365233 -1 68 8.132048 0.785914 1 69 8.245938 0.372366 1 70 6.543888 0.433164 1 71 -0.236713 -5.766721 -1 72 8.112593 0.295839 1 73 9.803425 1.495167 1 74 1.497407 -0.552916 -1 75 1.336267 -1.632889 -1 76 9.205805 -0.586480 1 77 1.966279 -1.840439 -1 78 8.398012 1.584918 1 79 7.239953 -1.764292 1 80 7.556201 0.241185 1 81 9.015509 0.345019 1 82 8.266085 -0.230977 1 83 8.545620 2.788799 1 84 9.295969 1.346332 1 85 2.404234 0.570278 -1 86 2.037772 0.021919 -1 87 1.727631 -0.453143 -1 88 1.979395 -0.050773 -1 89 8.092288 -1.372433 1 90 1.667645 0.239204 -1 91 9.854303 1.365116 1 92 7.921057 -1.327587 1 93 8.500757 1.492372 1 94 1.339746 -0.291183 -1 95 3.107511 0.758367 -1 96 2.609525 0.902979 -1 97 3.263585 1.367898 -1 98 2.912122 -0.202359 -1 99 1.731786 0.589096 -1 100 2.387003 1.573131 -1
5.2 SMO算法实现
定义数据结构体optStruct，用于缓存，提高运行速度。SMO算法具体实现如下（mySVMMLiA.py）
每个方法的作用，以及每行代码的作用，我都做了详细的注解，希望对大家的理解有帮助。
1 # -*- coding: utf-8 -*- 2 """ 3 Created on Tue Sep 4 16:58:16 2018 4 支持向量机代码实现 5 SMO(Sequential Minimal Optimization)最小序列优化 6 @author: weixw 7 """ 8 import numpy as np 9 #核转换函数（一个特征空间映射到另一个特征空间，低维空间映射到高维空间） 10 #高维空间解决线性问题，低维空间解决非线性问题 11 #线性内核 = 原始数据矩阵（100*2）与原始数据第一行矩阵转秩乘积（2*1） =>（100*1） 12 #非线性内核公式：k(x,y) = exp(-||x - y||**2/2*(e**2)) 13 #1.原始数据每一行与原始数据第一行作差， 14 #2.平方 15 def kernelTrans(dataMat, rowDataMat, kTup): 16 m,n=np.shape(dataMat) 17 #初始化核矩阵 m*1 18 K = np.mat(np.zeros((m,1))) 19 if kTup[0] == 'lin': #线性核 20 K = dataMat*rowDataMat.T 21 elif kTup[0] == 'rbf':#非线性核 22 for j in range(m): 23 #xi - xj 24 deltaRow = dataMat[j,:] - rowDataMat 25 K[j] = deltaRow*deltaRow.T 26 #1*m m*1 => 1*1 27 K = np.exp(K/(-2*kTup[1]**2)) 28 else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized') 29 return K 30 31 #定义数据结构体，用于缓存，提高运行速度 32 class optStruct: 33 def __init__(self, dataSet, labelSet, C, toler, kTup): 34 self.dataMat = np.mat(dataSet) #原始数据，转换成m*n矩阵 35 self.labelMat = np.mat(labelSet).T #标签数据 m*1矩阵 36 self.C = C #惩罚参数，C越大，容忍噪声度小，需要优化；反之，容忍噪声度高，不需要优化； 37 #所有的拉格朗日乘子都被限制在了以C为边长的矩形里 38 self.toler = toler #容忍度 39 self.m = np.shape(self.dataMat)[0] #原始数据行长度 40 self.alphas = np.mat(np.zeros((self.m,1))) # alpha系数，m*1矩阵 41 self.b = 0 #偏置 42 self.eCache = np.mat(np.zeros((self.m,2))) # 保存原始数据每行的预测值 43 self.K = np.mat(np.zeros((self.m,self.m))) # 核转换矩阵 m*m 44 for i in range(self.m): 45 self.K[:,i] = kernelTrans(self.dataMat, self.dataMat[i,:], kTup) 46 47 #计算原始数据第k项对应的预测误差 1*m m*1 =>1*1 48 #oS：结构数据 49 #k： 原始数据行索引 50 def calEk(oS, k): 51 #f(x) = w*x + b 52 fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) 53 Ek = fXk - float(oS.labelMat[k]) 54 return Ek 55 56 #在alpha有改变都要更新缓存 57 def updateEk(oS, k): 58 Ek = calEk(oS, k) 59 oS.eCache[k] = [1, Ek] 60 61 62 #第一次通过selectJrand()随机选取j,之后选取与i对应预测误差最大的j（步长最大） 63 def selectJ(i, oS, Ei): 64 #初始化 65 maxK = -1 #误差最大时对应索引 66 maxDeltaE = 0 #最大误差 67 Ej = 0 # j索引对应预测误差 68 #保存每一行的预测误差值 1相对于初始化为0的更改 69 oS.eCache[i] = [1,Ei] 70 #获取数据缓存结构中非0的索引列表(先将矩阵第0列转化为数组) 71 validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] 72 #遍历索引列表，寻找最大误差对应索引 73 if len(validEcacheList) > 1: 74 for k in validEcacheList: 75 if k == i: 76 continue 77 Ek = calEk(oS, k) 78 deltaE = abs(Ei - Ek) 79 if(deltaE > maxDeltaE): 80 maxK = k 81 maxDeltaE = deltaE 82 Ej = Ek 83 return maxK, Ej 84 else: 85 #随机选取一个不等于i的j 86 j = selectJrand(i, oS.m) 87 Ej = calEk(oS, j) 88 return j,Ej 89 90 #随机选取一个不等于i的索引 91 def selectJrand(i, m): 92 j = i 93 while (j == i): 94 j = int(np.random.uniform(0, m)) 95 return j 96 97 #alpha范围剪辑 98 def clipAlpha(aj, L, H): 99 if aj > H: 100 aj = H 101 if aj < L: 102 aj = L 103 return aj 104 105 #从文件获取特征数据，标签数据 106 def loadDataSet(fileName): 107 dataSet = []; labelSet = [] 108 fr = open(fileName) 109 for line in fr.readlines(): 110 #分割 111 lineArr = line.strip().split('\t') 112 dataSet.append([float(lineArr[0]), float(lineArr[1])]) 113 labelSet.append(float(lineArr[2])) 114 return dataSet, labelSet 115 116 #计算 w 权重系数 117 def calWs(alphas, dataSet, labelSet): 118 dataMat = np.mat(dataSet) 119 #1*100 => 100*1 120 labelMat = np.mat(labelSet).T 121 m, n = np.shape(dataMat) 122 w = np.zeros((n, 1)) 123 for i in range(m): 124 w += np.multiply(alphas[i]*labelMat[i], dataMat[i,:].T) 125 return w 126 #计算原始数据每一行alpha,b，保存到数据结构中，有变化及时更新 127 def innerL(i, oS): 128 #计算预测误差 129 Ei = calEk(oS, i) 130 #选择第一个alpha，违背KKT条件2 131 #正间隔，负间隔 132 if ((oS.labelMat[i] * Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.toler) and (oS.alphas[i] > 0)): 133 #第一次随机选取不等于i的数据项，其后根据误差最大选取数据项 134 j, Ej = selectJ(i, oS, Ei) 135 #初始化，开辟新的内存 136 alphaIold = oS.alphas[i].copy() 137 alphaJold = oS.alphas[j].copy() 138 #通过 a1y1 + a2y2 = 常量 139 # 0 <= a1,a2 <= C 求出L,H 140 if oS.labelMat[i] != oS.labelMat[j]: 141 L = max(0, oS.alphas[j] - oS.alphas[i]) 142 H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 143 else: 144 L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 145 H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 146 if L == H : 147 print ("L == H") 148 return 0 149 #内核分母 K11 + K22 - 2K12 150 eta = oS.K[i, i] + oS.K[j, j] - 2.0*oS.K[i, j] 151 if eta <= 0: 152 print ("eta <= 0") 153 return 0 154 #计算第一个alpha j 155 oS.alphas[j] += oS.labelMat[j]*(Ei - Ej)/eta 156 #修正alpha j的范围 157 oS.alphas[j] = clipAlpha(oS.alphas[j], L, H) 158 #alpha有改变，就需要更新缓存数据 159 updateEk(oS, j) 160 #如果优化后的alpha 与之前的alpha变化很小，则舍弃，并重新选择数据项的alpha 161 if (abs(oS.alphas[j] - alphaJold) < 0.00001): 162 print ("j not moving enough, abandon it.") 163 return 0 164 #计算alpha对的另一个alpha i 165 # ai_new*yi + aj_new*yj = 常量 166 # ai_old*yi + ai_old*yj = 常量 167 # 作差=> ai = ai_old + yi*yj*(aj_old - aj_new) 168 oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) 169 #alpha有改变，就需要更新缓存数据 170 updateEk(oS, i) 171 #计算b1,b2 172 # y(x) = w*x + b => b = y(x) - w*x 173 # w = aiyixi(i= 1->N求和) 174 #b1_new = y1_new - (a1_new*y1*k11 + a2_new*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量)) 175 #b1_old = y1_old - (a1_old*y1*k11 + a2_old*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量)) 176 #作差=> b1_new = b1_old + (y1_new - y1_old) - y1*k11*(a1_new - a1_old) - y2*k21*(a2_new - a2_old) 177 # => b1_new = b1_old + Ei - yi*(ai_new - ai_old)*kii - yj*(aj_new - aj_old)*kij 178 #同样可推得 b2_new = b2_old + Ej - yi*(ai_new - ai_old)*kij - yj*(aj_new - aj_old)*kjj 179 bi = oS.b - Ei - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[i,j] 180 bj = oS.b - Ej - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,j] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[j,j] 181 #首选alpha i，相对alpha j 更准确 182 if (0 < oS.alphas[i]) and (oS.alphas[i] < oS.C): 183 oS.b = bi 184 elif (0 < oS.alphas[j]) and (oS.alphas[j] < oS.C): 185 oS.b = bj 186 else: 187 oS.b = (bi + bj)/2.0 188 return 1 189 else: 190 return 0 191 192 #完整SMO核心算法，包含线性核核非线性核，返回alpha,b 193 #dataSet 原始特征数据 194 #labelSet 标签数据 195 #C 凸二次规划参数 196 #toler 容忍度 197 #maxInter 循环次数 198 #kTup 指定核方式 199 #程序逻辑： 200 #第一次全部遍历，遍历后根据alpha对是否有修改判断， 201 #如果alpha对没有修改，外循环终止；如果alpha对有修改，则继续遍历属于支持向量的数据。 202 #直至外循环次数达到maxIter 203 #相比简单SMO算法，运行速度更快，原因是： 204 #1.不是每一次都全量遍历原始数据，第一次遍历原始数据， 205 #如果alpha有优化，就遍历支持向量数据，直至alpha没有优化，然后再转全量遍历，这是如果alpha没有优化，循环结束； 206 #2.外循环不需要达到maxInter次数就终止； 207 def smoP(dataSet, labelSet, C, toler, maxInter, kTup = ('lin', 0)): 208 #初始化结构体类，获取实例 209 oS = optStruct(dataSet, labelSet, C, toler, kTup) 210 iter = 0 211 #全量遍历标志 212 entireSet = True 213 #alpha对是否优化标志 214 alphaPairsChanged = 0 215 #外循环 终止条件：1.达到最大次数 或者 2.alpha对没有优化 216 while (iter < maxInter) and ((alphaPairsChanged > 0) or (entireSet)): 217 alphaPairsChanged = 0 218 #全量遍历 ，遍历每一行数据 alpha对有修改，alphaPairsChanged累加 219 if entireSet: 220 for i in range(oS.m): 221 alphaPairsChanged += innerL(i, oS) 222 print ("fullSet, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged)) 223 iter += 1 224 else: 225 #获取(0，C)范围内数据索引列表，也就是只遍历属于支持向量的数据 226 nonBounds = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 227 for i in nonBounds: 228 alphaPairsChanged += innerL(i, oS) 229 print ("non-bound, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged)) 230 iter += 1 231 #全量遍历->支持向量遍历 232 if entireSet: 233 entireSet = False 234 #支持向量遍历->全量遍历 235 elif alphaPairsChanged == 0: 236 entireSet = True 237 print ("iteation number: %d"% iter) 238 print ("entireSet :%s"% entireSet) 239 print ("alphaPairsChanged :%d"% alphaPairsChanged) 240 return oS.b,oS.alphas 241 242 #绘制支持向量 243 def drawDataMap(dataArr,labelArr,b,alphas): 244 import matplotlib.pyplot as plt 245 #alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用 246 svInd=np.nonzero(alphas.A>0)[0] 247 #分类数据点 248 classified_pts = {'+1':[],'-1':[]} 249 for point,label in zip(dataArr,labelArr): 250 if label == 1.0: 251 classified_pts['+1'].append(point) 252 else: 253 classified_pts['-1'].append(point) 254 fig = plt.figure() 255 ax = fig.add_subplot(111) 256 #绘制数据点 257 for label,pts in classified_pts.items(): 258 pts = np.array(pts) 259 ax.scatter(pts[:, 0], pts[:, 1], label = label) 260 #绘制分割线 261 w = calWs(alphas, dataArr, labelArr) 262 #函数形式：max( x ,key=lambda a : b ) # x可以是任何数值，可以有多个x值 263 #先把x值带入lambda函数转换成b值，然后再将b值进行比较 264 x1, _=max(dataArr, key=lambda x:x[0]) 265 x2, _=min(dataArr, key=lambda x:x[0]) 266 a1, a2 = w 267 y1, y2 = (-b - a1*x1)/a2, (-b - a1*x2)/a2 268 #矩阵转化为数组.A 269 ax.plot([x1, x2],[y1.A[0][0], y2.A[0][0]]) 270 271 #绘制支持向量 272 for i in svInd: 273 x, y= dataArr[i] 274 ax.scatter([x], [y], s=150, c ='none', alpha=0.7, linewidth=1.5, edgecolor = '#AB3319') 275 plt.show() 276 277 #alpha>0对应的数据才是支持向量，过滤不是支持向量的数据 278 sVs= np.mat(dataArr)[svInd] #get matrix of only support vectors 279 print ("there are %d Support Vectors.\n" % np.shape(sVs)[0]) 280 281 #训练结果 282 def getTrainingDataResult(dataSet, labelSet, b, alphas, k1=1.3): 283 datMat = np.mat(dataSet) 284 #100*1 285 labelMat = np.mat(labelSet).T 286 #alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用 287 svInd=np.nonzero(alphas.A>0)[0] 288 sVs=datMat[svInd] 289 labelSV = labelMat[svInd]; 290 m,n = np.shape(datMat) 291 errorCount = 0 292 for i in range(m): 293 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 294 # y(x) = w*x + b => b = y(x) - w*x 295 # w = aiyixi(i= 1->N求和) 296 predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b 297 if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1 298 print ("the training error rate is: %f" % (float(errorCount)/m)) 299 300 def getTestDataResult(dataSet, labelSet, b, alphas, k1=1.3): 301 datMat = np.mat(dataSet) 302 #100*1 303 labelMat = np.mat(labelSet).T 304 #alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用 305 svInd=np.nonzero(alphas.A>0)[0] 306 sVs=datMat[svInd] 307 labelSV = labelMat[svInd]; 308 m,n = np.shape(datMat) 309 errorCount = 0 310 for i in range(m): 311 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 312 # y(x) = w*x + b => b = y(x) - w*x 313 # w = aiyixi(i= 1->N求和) 314 predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b 315 if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1 316 print ("the test error rate is: %f" % (float(errorCount)/m)) 317 318
5.3 测试代码（testMySVMMLiA.py）
1 # -*- coding: utf-8 -*- 2 """ 3 Created on Wed Sep 5 15:22:26 2018 4 5 @author: weixw 6 """ 7 8 9 import mySVMMLiA as sm 10 11 #通过训练数据计算 b, alphas 12 dataArr,labelArr = sm.loadDataSet('trainingData.txt') 13 b, alphas = sm.smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 0.10)) 14 sm.drawDataMap(dataArr,labelArr,b,alphas) 15 sm.getTrainingDataResult(dataArr, labelArr, b, alphas, 0.10) 16 dataArr1,labelArr1 = sm.loadDataSet('testData.txt') 17 #测试结果 18 sm.getTestDataResult(dataArr1, labelArr1, b, alphas, 0.10)
通过训练数据计算出 b， 权重矩阵，从而分类超平面和决策分类函数就明确了，然后测试数据以决策分类函数进行预测。
这里采用高斯核RBF。
5.4 运行结果
“j not moving enough, abandon it” 表示数据项对应的 和
非常接近，不需要优化；
“fullSet” 表示全量数据遍历；
“non-bound” 表示非边界遍历，也就是只遍历属于支持向量的数据项。
另外我将支持向量的数据项绘制出来了，这样更直观。
可以看出，有77个支持向量，训练差错率是0，测试差错率6%。
6. 参考文献
《统计学习方法》第 七 章
《机器学习实战》第 六 章
《从零推导支持向量机》
博客：https://www.cnblogs.com/bentuwuying/p/6444516.html
博客：http://www.cnblogs.com/jerrylead/archive/2011/03/18/1988419.html
博客：https://www.cnblogs.com/zhizhan/p/4430253.html
博客：https://www.cnblogs.com/xxrxxr/p/7538430.html
不要让懒惰占据你的大脑，不要让妥协拖垮了你的人生。青春就是一张票，能不能赶上时代的快车，你的步伐就掌握在你的脚下。
