本文系作者原创,转载请注明出处:https://cloud.tencent.com/developer/article/1691209
最基本的应用是数据分类,特别是对于非线性不可分数据集。支持向量机不仅能对非线性可分数据集进行分类,对于非线性不可分数据集的也可以分类
(我认为这才是支持向量机的真正魅力所在,因为现实场景中,样本数据往往是非线性不可分的)。
现实场景一 :样本数据大部分是线性可分的,但是只是在样本中含有少量噪声或特异点,去掉这些噪声或特异点后线性可分 => 用支持向量机的软间隔方法进行分类;
现实场景二:样本数据完全线性不可分 => 引入核函数,将低维不可分的非线性数据集转化为高维可分的数据集,用支持向量机的软间隔方法进行分类;
一定要明白一点:对分类器准确性有影响只是样本中的支持向量,因此,其他样本在计算分类器的权重矩阵时可以直接过滤掉,大大节省运行时间。
目标:获取离边界最近样本点到超平面最远距离
博客 https://www.cnblogs.com/zhizhan/p/4430253.html 关于SVM的演化说得很透彻,也很形象,这里借用一下。
SVM中最关键的思想之一就是引入和定义了“间隔”这个概念。这个概念本身很简单,以二维空间为例,就是点到分类直线之间的距离。假设直线为y=wx+b,那么只要使所有正分类点到该直线的距离与所有负分类点到该直线的距离的总和达到最大,这条直线就是最优分类直线。这样,原问题就转化为一个约束优化问题,可以直接求解。这叫做硬间隔最大化,得到的SVM模型称作硬间隔支持向量机。
但是新问题出现了,在实际应用中,我们得到的数据并不总是完美的线性可分的,其中可能会有个别噪声点,他们错误的被分类到了其他类中。如果将这些特异的噪点去除后,可以很容易的线性可分。但是,我们对于数据集中哪些是噪声点却是不知道的,如果以之前的方法进行求解,会无法进行线性分开。是不是就没办法了呢?假设在y=x+1直线上下分为两类,若两类中各有对方的几个噪点,在人的眼中,仍然是可以将两类分开的。这是因为在人脑中是可以容忍一定的误差的,仍然使用y=x+1直线分类,可以在最小误差的情况下进行最优的分类。同样的道理,我们在SVM中引入误差的概念,将其称作“松弛变量”。通过加入松弛变量,在原距离函数中需要加入新的松弛变量带来的误差,这样,最终的优化目标函数变成了两个部分组成:距离函数和松弛变量误差。这两个部分的重要程度并不是相等的,而是需要依据具体问题而定的,因此,我们加入权重参数C,将其与目标函数中的松弛变量误差相乘,这样,就可以通过调整C来对二者的系数进行调和。如果我们能够容忍噪声,那就把C调小,让他的权重降下来,从而变得不重要;反之,我们需要很严格的噪声小的模型,则将C调大一点,权重提升上去,变得更加重要。通过对参数C的调整,可以对模型进行控制。这叫做软间隔最大化,得到的SVM称作软间隔支持向量机。
之前的硬间隔支持向量机和软间隔支持向量机都是解决线性可分数据集或近似线性可分数据集的问题的。但是如果噪点很多,甚至会造成数据变成了线性不可分的,那该怎么办?最常见的例子是在二维平面笛卡尔坐标系下,以原点(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是在对问题不了解的情况下,对最广泛问题效果都很不错的核函数。因此,使用范围也最广。
求使几何间隔最大的分离超平面
以及相应的分类决策函数
也就是求出w,b
为了求解线性可分支持向量机的最优化问题,将它作为原始最优化问题,应用拉格朗日对偶性,通过求解对偶问题得到原始问题最优解。
这就是线性可分支持向量机的对偶算法,这样做的优点,一是对偶问题往往更容易求解;二是自然引入核函数,进而推广到非线性分类问题。
的求解就转化为
的求解。上述公式就是对不同变量求偏导,具体推导过程见《统计学习方法》第7章 或者相关博客。
上图的总结有点小问题,
分类不能合并,正确的分类描述应该是:
在超平面上或误分的样本点是不能正确分类的。
高效实现支持向量机的算法是SMO算法,其基本思路是:如果所有变量的解都满足此最优化的KKT条件,那么这个最优化问题的解就得到了。因为KKT条件是该优化问题的充分必要条件。否则,
选择两个变量,固定其他变量,针对这两个变量构建一个二次规划问题,这个二次规划问题关于这两个变量的解应该更接近原始二次规划问题的解,因为这会使原始二次规划问题的目标函数值变得更小。
重要的是,这时子问题可以通过解析方法求解,这样可以大大提高整个算法的计算速度。
子问题有两个变量,一个是违反KKT条件最严重的那一个,另一个由约束条件自动确定。如此,SMO算法将原问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
由
并且
推出
范围
由
假定
不固定,其他拉格朗日因子固定(也就是常量),得到
常量 ,
为变量,
为因变量,L,H的范围就是
的 范围
,
第一个
公式如下,推导过程见《统计学习方法》第7章或者相关博客。
选择一个违背KKT的数据项,根据下面结论得到优化后
,
s表示
,从而可以得到另一个优化后的
,同样也需要进行约束。
在每次完成两个变量的优化后,都要重新计算阈值b。具体推导见《统计学习方法》
这里使用两个原始数据文件 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
定义数据结构体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
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。
“j not moving enough, abandon it” 表示数据项对应的
和
非常接近,不需要优化;
“fullSet” 表示全量数据遍历;
“non-bound” 表示非边界遍历,也就是只遍历属于支持向量的数据项。
另外我将支持向量的数据项绘制出来了,这样更直观。
可以看出,有77个支持向量,训练差错率是0,测试差错率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
不要让懒惰占据你的大脑,不要让妥协拖垮了你的人生。青春就是一张票,能不能赶上时代的快车,你的步伐就掌握在你的脚下。