机器学习-支持向量机(python3代码实现)
支持向量机
哈尔滨工程大学-537
算法原理：
一、寻找最大间隔
如下图所示，用一条分割线将两类点分割开来(二维的是一条分割线，多维的就是分隔面)，显然三条线都能将两类点分割开来，然而，从直观来看，红色的分割线显然分割效果最好。为什么这么说呢？
因为红色的分割线到两边最近的点的距离更远。可以直观把两边的两类点想象成地雷，我们有一支红军要通过这片雷区，显然，沿着绿色和灰色的路线行军，两边不会踩到地雷的安全区域非常的窄，而沿着红色的路线行军，安全区域明显更宽。我们的目的就是要找到能够使安全区域最宽的行军路线，最终使红军夺得革命战争的胜利。
如下图所示，距离分割线最近的几个点的位置，决定了分割线的位置，形象的来说，就是距离红军队伍最近的几个地雷的位置，决定了红军穿过雷区的行军路线，如果距离红军最近的地雷的位置发生改变，红军的行军路线就必须随之改变，否则安全区域就有可能变窄，踩到地雷的可能性就会增加。
那么此时要解决的目标就非常明确了，即找到距离分割线最近的那几个点，而这几个距离分割线最近的点，就叫做支持向量，这就是支持向量机一词的由来。
那么如何找到距离分割线最近的那几个点？
由高中数学可知：空间中一个点(x1,y1,z1)(x1,y1,z1)(x_1,y_1,z_1)到平面线ax+by+cz+d=0ax+by+cz+d=0ax+by+cz+d=0 的距离为：ax1+by1+cz1+da2+b2+c2√ax1+by1+cz1+da2+b2+c2\frac{ax_1+by_1+cz_1+d}{\sqrt{a^2+b^2+c^2}} ；
扩展为多维的情况，点Xi=(xi1,xi2...xik)Xi=(xi1,xi2...xik)X^i=(x^i{_1}, x^i{_2}...x^i{_k})到一个超平面WTX+b=0WTX+b=0W{^T}X+b=0的距离为：WTXi+b||W||WTXi+b||W||\frac{W{^T}X{^i}+b}{||W||}，其中WWW和X" role="presentation" style="position: relative;">XXX为kkk维向量(因为Xi" role="presentation" style="position: relative;">XiXiX^i有kkk个特征)。
于是当前的任务就是要找到WTXi+b||W||" role="presentation" style="position: relative;">WTXi+b||W||WTXi+b||W||\frac{W{^T}X{^i}+b}{||W||}值最小的数据点，将该点的WTXi+b||W||WTXi+b||W||\frac{W{^T}X{^i}+b}{||W||}最大化，此时的WWW和b" role="presentation" style="position: relative;">bbb就是我们要找的最优分割超平面的参数。
由高中数学可知：若点(x1,y1)(x1,y1)(x_1,y_1)在直线y=ax+by=ax+by=ax+b的上侧，则将点(x1,y1)(x1,y1)(x_1,y_1)带入直线得ax1+b−y1>0ax1+b−y1>0ax_1+b-y_1>0，反之，若在下侧，则带入直线得ax1+b−y1<0ax1+b−y1<0ax_1+b-y_1;
推广到多维的情况：若数据点XiXiX^i在超平面正侧，WTXi+b>0WTXi+b>0{W{^T}X{^i}+b}>0，那么将在这一侧的数据点定义为1类，即类别标签yiyiy^i为1，那么yi(WTXi+b)>0yi(WTXi+b)>0y{^i}{(W{^T}X{^i}+b)}>0；
反之，若数据点XiXiX^i在超平面的负侧，WTXi+b<0WTXi+b<0{W{^T}X{^i}+b},那么将这一侧的数据点定义为-1类，即类别标签yiyiy^i为-1，那么yi(WTXi+b)>0yi(WTXi+b)>0y{^i}{(W{^T}X{^i}+b)}>0 ， 这样在比较大小的时候，就避免了负数的出现。
那么此时，首要任务就是找到yi(WTXi+b)||W||yi(WTXi+b)||W||\frac{y^i(W{^T}X{^i}+b)}{||W||}最小的数据点，并将该点的yi(WTXi+b)||W||yi(WTXi+b)||W||\frac{y^i(W{^T}X{^i}+b)}{||W||}值最大化。
若限制yi(WTXi+b)≥1yi(WTXi+b)≥1y^i(W{^T}X{^i}+b)\geq1，则距离超平面最近的点的yi(WTXi+b)yi(WTXi+b)y^i(W{^T}X{^i}+b)应等于1，而||W||||W||||W||则越大，说明该点离超平面越近。
如下图，可以更加直观的理解以上说法，虚线WTX+b−1=0WTX+b−1=0W^TX+b-1=0和虚线WTX+b+1=0WTX+b+1=0W^TX+b+1=0分别是两条与直线WTX+b=0WTX+b=0W^TX+b=0平行的直线(由高中数学可知,WTX+b−1=0WTX+b−1=0W^TX+b-1=0在直线上侧，WTX+b+1=0WTX+b+1=0W^TX+b+1=0在直线下侧)，通过归一化系数W，可以使最后的常数一直保持+1和-1，也就是说，这两条虚线可以在平面上任意移动，而始终保持WTX+b−1=0WTX+b−1=0W^TX+b-1=0和WTX+b+1=0WTX+b+1=0W^TX+b+1=0的形式。那么现在的要求就是，让这两条虚线之间的距离最大，且要保证所有点都在这两条虚线之外(或在虚线之上)，即1类样本点都在WTX+b−1=0WTX+b−1=0W^TX+b-1=0的正侧(或在线上)，即WTXi+b−1≥0WTXi+b−1≥0W^TX_i+b-1\geq0；而-1类样本点都在WTX+b+1=0WTX+b+1=0W^TX+b+1=0的负侧(或在线上)，即WTXi+b+1≤0WTXi+b+1≤0W^TX_i+b+1\leq0。
如果把+1和-1项都移到等式右边，就分别为：WTXi+b≥+1WTXi+b≥+1W^TX_i+b\geq+1和WTXi+b≤−1WTXi+b≤−1W^TX_i+b\leq-1，再分别乘上它们的类别标签yiyiy^i，就可以用一个式子来表示了，即WTXi+b≥1WTXi+b≥1W^TX_i+b\geq1。
那么此时的首要任务可以这样表述：给定训练样本(Xi,yi)(Xi,yi)(X^i,y^i)，在yi(WTXi+b)≥1yi(WTXi+b)≥1y^i(W{^T}X{^i}+b)\geq1的限制条件下，使yi(WTXi+b)||W||yi(WTXi+b)||W||\frac{y^i(W{^T}X{^i}+b)}{||W||}尽可能的大。
而能对分割超平面产生影响的只有支持向量，即位于虚线上的样本点，即yi(WTXi+b)=1yi(WTXi+b)=1y^i(W{^T}X{^i}+b)=1的样本点，那么此时只需要使分母，即||W||||W||||W||尽可能的大，为了后面计算方便，改写为使12||W||12||W||\frac12||W||即12WTW12WTW\frac12W^TW尽可能的大。此时的WWW和b" role="presentation" style="position: relative;">bbb就是最优超平面的参数。
由高数内容可知：带限制条件的求极值，可以应用拉格朗日乘子法。即求12WTW12WTW\frac12W^TW的极值，限制条件为1−yi(WTXi+b)≤01−yi(WTXi+b)≤01-y^i(W{^T}X{^i}+b)\leq0(拉格朗日乘子法要求限制条件为f(x)≤0f(x)≤0f(x)\leq0的形式)。由于对于每一个数据点Xi(i=1,2...n)Xi(i=1,2...n)X^i(i=1,2...n)，都作此限定条件，即任意的1−yi(WTXi+b)≤01−yi(WTXi+b)≤01-y^i(W{^T}X{^i}+b)\leq0，(i=1,2...n)(i=1,2...n)(i=1,2...n)。
于是由拉格朗日乘子法得：
L(W,\alpha,,b)=\frac{1}{2}W^TW+\sum_{i=1}^n\alpha_i(1-y^i(W^TX^i+b))
要求什么样的WWW和b" role="presentation" style="position: relative;">bbb能使L(W,α,b)L(W,α,b)L(W,\alpha,b)取得极小值，应分别对WWW和b" role="presentation" style="position: relative;">bbb求导，使导数为0，于是得到：
\frac{\partial{L}}{\partial{W}}=W-\sum_{i=1}^n\alpha_iy^iX^i=0
\frac{\partial{l}}{\partial{b}}=\sum_{i=1}^n\alpha_iy^i=0
于是可得：
W=\sum_{i=1}^n\alpha_iy^iX^i
首先将L(W,α,b)L(W,α,b)L(W,\alpha,b)逐项展开，得到：
L(W,\alpha,b)=\frac12W^TW+\sum_{i=1}^n\alpha_i-\sum_{i=1}^n\alpha_iy^iW^TX_i-b\sum_{i=1}^n\alpha_iy^i
再将 WWW带回到L(W,α,b)" role="presentation" style="position: relative;">L(W,α,b)L(W,α,b)L(W,\alpha,b)中，并且 ∑ni=1αiyi=0∑i=1nαiyi=0\sum_{i=1}^n\alpha_iy^i=0，得到：
\frac12(\sum_{i=1}^n\alpha_iy^iX^i)^T(\sum_{i=1}^n\alpha_iy^iX^i)+\sum_{i=1}^n\alpha_i-\sum_{i=1}^n\alpha_iy^i(\sum_{i=1}^n\alpha_iy^iX^i)^TX_i
整理后可得：
J(\alpha)=\sum_{i=1}^n\alpha_i-\frac12\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jy^iy^j(X^i)^TX^i
于是得到此时的目标函数为J(α)J(α)J(\alpha)。即这样的WWW和b" role="presentation" style="position: relative;">bbb能够使得L(W,α,b)L(W,α,b)L(W,\alpha,b)取得极小值，并且将WWW和b" role="presentation" style="position: relative;">bbb带入L(W,α,b)L(W,α,b)L(W,\alpha,b)后得到J(α)J(α)J(\alpha)，现在的目标是将J(α)J(α)J(\alpha)最大化。并且限制条件为：
\sum_{i=1}^n\alpha_iy^i=0
\alpha_i\geq0
式子 J(α)J(α)J(\alpha)中的 Xi,yiXi,yiX^i,y^i是全体样本点，然而 J(α)J(α)J(\alpha)取得极大值仅仅依赖距离超平面最近的点，即只和支持向量有关。所以 J(α)J(α)J(\alpha)的解必然是稀疏的，即对于支持向量， α>0α>0\alpha>0，而对于其他样本点， α=0α=0\alpha=0。
故将样本点带入J(α)J(α)J(\alpha)，再分别对αiαi\alpha_i求导，使导数等于0，再根据∑ni=1αiyi=0∑i=1nαiyi=0\sum_{i=1}^n\alpha_iy^i=0,求解αiαi\alpha_i，且αiαi\alpha_i需满足限定条件αi≥0αi≥0\alpha_i\geq0。
若解不满足限定条件，则解可能在限制边界上，则将将其中一个αα\alpha以0带入，求解其他αα\alpha值
再根据式子：W=∑ni=1αiyiXiW=∑i=1nαiyiXiW=\sum_{i=1}^n\alpha_iy^iX^i，求解出WWW。
最后根据式子yi(WTXi+b)=1" role="presentation" style="position: relative;">yi(WTXi+b)=1yi(WTXi+b)=1y^i(W{^T}X{^i}+b)=1，求解出bbb。
于是便得到了最优超平面的参数W" role="presentation" style="position: relative;">WWW和bbb。
二、引入松弛变量
以上是针对比较完美的情况，即两类点100%线性可分，但是如果数据点并没有那么“干净”，这时候就可以引入所谓“松弛变量”，来允许有些数据点可以处于分割面的错误一侧，此时，限制条件为：
y^i(W^TX^i+b)\geq1-\xi_i
\xi_i\geq0
之前是使得 12WTW12WTW\frac12W^TW取极小值，现在是是使 12WTW+C∑ni=1ξi12WTW+C∑i=1nξi\frac12W^TW+C\sum_{i=1}^n\xi_i取得极小值。这里的 CCC是可调节的参数，当C" role="presentation" style="position: relative;">CCC比较大时，显然 ξiξi\xi_i比较小； CCC比较小时，则ξi" role="presentation" style="position: relative;">ξiξi\xi_i比较大。
于是此时应用拉格朗日乘子法可得：
L(W,b,\xi_i,\alpha,\mu)=\frac12W^TW+C\sum_{i=1}^n\xi_i+\sum_{i=1}^n\alpha_i(1-\xi_i-y^i(W^TX^i+b))-\sum_{i=1}^n\mu_i\xi_i
(这里将 ξi≥0ξi≥0\xi_i\geq0变成 −ξi≤0−ξi≤0-\xi_i\leq0的形式，以适用拉格朗日乘子法)
与之前同样的做法，求L(W,b,ξi,α,μ)L(W,b,ξi,α,μ)L(W,b,\xi_i,\alpha,\mu)取得极小值，应对WWW，b" role="presentation" style="position: relative;">bbb，以及引入的“松弛变量”ξiξi\xi_i分别求偏导(i=1,2...n)(i=1,2...n)(i=1,2...n)，可以看出对WWW和b" role="presentation" style="position: relative;">bbb求偏导结果没有改变，因此仍然可得：
W=\sum_{i=1}^n\alpha_iy^iX^i以及
\sum_{i=1}^n\alpha_iy^i=0带入整理后仍然可得：
J(\alpha)=\sum_{i=1}^n\alpha_i-\frac12\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jy^iy^j(X^i)^TX^i
而对ξi(i=1,2...n)ξi(i=1,2...n)\xi_i(i=1,2...n)分别求偏导得到：
C-\alpha_i-\mu_i=0
此时的任务仍然是使 J(α)J(α)J(\alpha)最大化，而此时的限制条件为：
\sum_{i=1}^n\alpha_i=0
\alpha_i\geq0
\mu_i\geq0
C-\alpha_i-\mu_i=0
后三个限制条件可以合并为一个限制条件：
0\leq\alpha_i\leq C
于是此时便可以与之前同样的解法求解 αiαi\alpha_i，唯一不同的是， αiαi\alpha_i的值被限定为小于一个设定的常数C。
三、SMO求解算法(简化版)
现在的问题是如何求解αi(i=1,2...n)αi(i=1,2...n)\alpha_i(i=1,2...n)，(显然有多少个样本点，就有多少个αα\alpha)，大神Platt给出了SMO算法进行求解，这里可以跟随实际的python代码，一行一行的分析SMO算法的原理，由于代码较复杂，在代码中敲大量注释影响代码清晰性，所以在代码之外进行相应的解释。
1）该段代码进行需要的各种库的导入，最后一行%matplotlib inline是为了保证在jupyter notebook 的运行环境下能够在浏览器上输出绘制的图像。
import numpy as np
import random
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline 2）然后用pandas读取txt文件，并调用head方法显示文件的前几行。
file_name = r"D:\python_data\MachineLearningInaction\machinelearninginaction\Ch06\testSet.txt"
file = pd.read_table(file_name,header=None, names=["factor1","factor2","class"])
file.head()3）可以看到所读取的txt文件的前几行如下图，每个样本有2个特征，类标签分别为+1或-1。
4）接下来再绘制一张样本点分布的散点图，蓝色的是正例，即+1类的点；红色的是负例，即-1类的点。横坐标为“factor1”，纵坐标为“factor2”。
positive = file[file["class"] == 1]
negative = file[file["class"] == -1]
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(positive["factor1"], positive["factor2"], s=30, c="b", marker="o", label="class 1")
ax.scatter(negative["factor1"], negative["factor2"], s=30, c="r", marker="x", label="class -1")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")5）运行后，在浏览器上输出了如下的图像。
6）现在对数据已经有了清晰直观的了解，接下来需要把pandas读取出来的整个表格进行切分，前2列作为样本点的特征矩阵，最后一列作为类标签向量。下面函数输入的是pandas都取出来的文件file，调用as_matrix方法将其转化为矩阵orig_data，将这个矩阵去除最后一列的部分作为样本点的特征矩阵data_mat，而最后一列作为类标签向量label_mat。
def load_data_set(file):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_mat7）接下来调用该函数，得到样本点的特征矩阵data_mat和对应的类标签向量label_mat。
data_mat, label_mat = load_data_set(file) 8）还需要定义一个select_jrand函数，该函数的功能是给定输入i和m，随机输出一个0到m之间的与i不同的整数。这个函数在后面将会用到。
def select_jrand(i, m):j = iwhile(j==i):j = int(random.uniform(0,m))return j9）此外还需要定义一个clip_alpha函数，该函数的功能是输入将aj限制在L和H之间
def clip_alpha(aj, H, L):if aj > H:aj = Hif L > aj:aj = Lreturn aj10）接下来是主体功能的实现函数，由于这个函数太长，为了方便理解，将以下代码分成一小段一小段的，这样就可以理解每一小段的功能。
10.1）循环外的初始化工作：首先该函数的输入为特征矩阵、类标签向量，常数C，容错率，最大迭代次数。首先将输入的特征矩阵和类标签向量都转化为矩阵(若特征矩阵和类标签向量已经是矩阵，则不需要这个步骤，若输入的是列表的形式，则需要该步骤)，设置b为0，m和n分别为特征矩阵的行数和列数。创建一个m行1列的αα\alpha向量(即有多少个样本点就有多少个αα\alpha)，迭代次数初始化为0(每次迭代在αα\alpha向量中寻找一对可以改变的αα\alpha，若没有找到则迭代次数+1，达到最大迭代次数还没有找到可以改变的αα\alpha，说明已经αα\alpha已经优化的差不多了，则退出循环)。
10.2）内循环的初始化工作：现在开始第一轮迭代，将改变的一对αα\alpha的数量alpha_pairs_changed初始化为0，然后遍历m个αα\alpha依次作为第一个αα\alpha，再在(0,m)范围内随机的找到第二个αα\alpha(第二个αα\alpha与第一个不能是同一个，这就用到了之前定义的select_jrand函数)。
10.3)第一小段代码：首先从i等于0开始，已知公式W=∑niαiyiXiW=∑inαiyiXiW=\sum_i^n\alpha_iy^iX^i，若第0号数据点是支持向量，则该数据点必然在虚线WTX+b+1=0WTX+b+1=0W^TX+b+1=0或WTX+b−1=0WTX+b−1=0W^TX+b-1=0上，即WTX+b=−1WTX+b=−1W^TX+b=-1或WTX+b=+1WTX+b=+1W^TX+b=+1。
先将0号数据点带入公式求出WWW，即代码中的WT_i，再将0号数据点和求得的W" role="presentation" style="position: relative;">WWW带入WTX+bWTX+bW^TX+b求出WTX+bWTX+bW^TX+b的值，即代码中的f_xi。
代码中的E_i为得到的f_xi与实际的WTX+bWTX+bW^TX+b之间的误差，而这个误差大有说道，最完美的情况是若数据点是支持向量，则将数据点带入WTX+bWTX+bW^TX+b应该等于-1或+1，即该数据点应该在虚线WTX+b+1=0WTX+b+1=0W^TX+b+1=0或虚线WTX+b−1=0WTX+b−1=0W^TX+b-1=0上，而如果该数据点的类标签yiyiy^i和误差E_i的乘积大于容错率，说明该数据点属于+1类且在WTX+b−1=0WTX+b−1=0W^TX+b-1=0的上侧，或者该数据点属于-1类且在WTX+b+1=0WTX+b+1=0W^TX+b+1=0的下侧。也就是说该点在属于自己的那一类的群体里，但是太靠后了且靠后的超过限度了(即大于容错率)，根本无法作为支持向量，而这个点对应的αα\alpha值又大于0，说明这明明就不该是靠后站的，那么就要调整这个数据点的αα\alpha了。而另一种情况是该数据点的类标签yiyiy^i和误差E_i的乘积小于负容错率，说明该点太靠前了已经处于虚线WTX+b+1WTX+b+1W^TX+b+1和虚线WTX+b−1WTX+b−1W^TX+b-1之间了，且该点的αα\alpha值又小于常数C，说明该点也不应该是靠前的，那就要调整这个数据点的αα\alpha了
若经过上述判断，决定要调整该数据点的αα\alpha，那么就调用之前定义的select_jrand函数随机的选择另一个要调整的αα\alpha，对另一个alphaalphaalpha对应的数据点也计算相应的误差E_j，为了记录αα\alpha的调整情况，分别copy一下作为αiαi\alpha_i和αjαj\alpha_j的旧值。
10.4）第二小段代码：由于有一个限制条件即∑niαiyi=0∑inαiyi=0\sum_i^n\alpha_iy^i=0，且任意的αα\alpha在(0,C)之间，所以当调整其中一个αiαi\alpha_i的时候，另一个αjαj\alpha_j的改变必须有一定限制，即yiαi+yjαj=kyiαi+yjαj=ky^i\alpha_i+y^j\alpha_j=k。即如下图：αiαi\alpha_i和αjαj\alpha_j必须在两条线段上。
根据上图，可以得到αjαj\alpha_j的取值范围(L,H)。
10.5）第三小段代码：下面为αjαj\alpha_j的更新公式：
\alpha_j = \alpha_j-\frac{y^j(E_i-E_j)}{\eta}
其中 η=2<Xi,Xj>−<Xi,Xi>−<Xj,Xj>η=2<Xi,Xj>−<Xi,Xi>−<Xj,Xj>\eta=2--。(尖括号表示两个向量的内积)
更新前先对 ηη\eta进行判断，若 η≥0η≥0\eta\geq0则退出本次循环，不进行更新。否则对 αjαj\alpha_j进行更新。更新后需要调用之前定义的clip_alpha函数，将更新后的 αjαj\alpha_j缩小在这个范围之内，最后检验一下 αjαj\alpha_j值改变了多少，如果改变的太少，就不要更新了，直接退出本次循环寻找下一对 αα\alpha，改变那点值没什么意义。
最后既然 αjαj\alpha_j更新了， αiαi\alpha_i也得随之改变，且该变量应该是相等的，根据之前所述的 yiαi+yjαj=kyiαi+yjαj=ky^i\alpha_i+y^j\alpha_j=k，可以得到更新后的 αiαi\alpha_i的值。
10.6）第四小段代码：对于样本点i，已经更新了其αiαi\alpha_i值，若b值也是正确的，则该点应该从更新之前的不在WTX+b±1=0WTX+b±1=0W^TX+b\pm1=0上，更新到位于WTX+b±1=0WTX+b±1=0W^TX+b\pm1=0上，亦即将样本点i带入，得到Enewi=WnewTXi+bnew−yi=0Einew=WnewTXi+bnew−yi=0E_i^{new}={W_{new}}^TX^i+b^{new}-{y^i}=0。
同时又已知Eoldi=WoldTXi+bold−yiEiold=WoldTXi+bold−yiE_i^{old}={W_{old}}^TX^i+b^{old}-y^i，根据上述两个公式，可以得到：
b^{new}=b^{old}-E_i^{old}+{W_{old}}^TX^i-{W_{new}}^TX^i
前两项不用多说，最后两项的相减，有些说道，因为W=∑ni=1αiyiXiW=∑i=1nαiyiXiW=\sum_{i=1}^n\alpha_iy^iX^i，所以其实WoldWoldW_{old}和WnewWnewW_{new}的区别只在里面的αiyiXiαiyiXi\alpha_iy^iX^i和αjyjXjαjyjXj\alpha_jy^jX^j不同，只要把公式展开成向量的形式写出来就很容易看出来，二者相减之后的结果就是:
y^i(\alpha{_i}^{old}-\alpha{_i}^{new})+y^j(\alpha_j^{old}-\alpha_j^{new})
将后两项带入，就可以得到 bnewbnewb^{new}。
当然这里得到的是根据样本点i求得的 bnewbnewb^{new}，根据样本点j做同样操作，可以得到另一个 bnewbnewb^{new}，此时若更新后的 αiαi\alpha_i在限定范围内，则 b=bnew1b=b1newb=b_1^{new},若 αjαj\alpha_j在限定范围内，则 b=bnew2b=b2newb=b_2^{new}，若 αiαi\alpha_i和 αjαj\alpha_j都在限定范围内，则 b=b1+b22b=b1+b22b=\frac{b_1+b_2}{2}。
10.7)收尾工作：若程序运行到此处，则已经更新了一对αα\alpha，并且也更新了bbb值，则将alpha_pairs_changed加上1，代表更新了一对，继续更新下一对。若程序没有执行更新α" role="presentation" style="position: relative;">αα\alpha的操作，即alpha_pairs_changed=0，则将iter+1，代表迭代了一次，否则将iter置为1，重新开始迭代。若程序迭代的次数达到了最大迭代次数，还没有αα\alpha被更新，说明已经更新的很好了，则可以退出while循环，返回b值和αα\alpha向量。
def smo_simple(data_mat, class_label, C, toler, max_iter):# 循环外的初始化工作data_mat = np.mat(data_mat); label_mat = np.mat(class_label) b = 0m,n = np.shape(data_mat) alphas = np.zeros((m,1)) iter = 0 while iter < max_iter:# 内循环的初始化工作alpha_pairs_changed = 0for 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 < -toler) and (alphas[i] < C)) or \((label_mat[i]*Ei > toler) and (alphas[i] > 0)): j = 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,:].Tif eta >= 0: continuealphas[j] = (alphas[j] - label_mat[j]*(Ei - Ej))/etaalphas[j] = clip_alpha(alphas[j], H, L) if (abs(alphas[j] - alpha_jold) < 0.00001):continuealphas[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 = b1elif (0 < alphas[j]) and (C > alphas[j]):b = b2else:b = (b1 + b2)/2.0alpha_pairs_changed += 1if (alpha_pairs_changed == 0): iter += 1else: iter = 0return b, alphas11）接下来调用该函数，得到b值和αα\alpha向量，由于αα\alpha的值大多为0，所以只打印大于零的αα\alpha值。
b, alphas = smo_simple(data_mat, label_mat, 0.6, 0.001, 10)
print(b, alphas[alphas>0])
12）辛辛苦苦走到这里，就看到几个所谓的αα\alpha和一个b值，真的不够。“锦瑟无端五十弦，一弦一柱思华年”，每一个算法的推导，每一段代码的理解，都给我留下了深刻的印象。过程的享受是一种润物无声的细腻，但我想要的更是成功之后那会“当凌绝顶，一览众山小”的快意，更要那“冲天香阵透长安，满城尽戴黄金甲”的震撼。于是，我又敲下了如下代码，得到了，那虽不太完美，却凝聚了我心血的分、割、线。
support_x = []
support_y = []
class1_x = []
class1_y = []
class01_x = []
class01_y = []
for i in range(100):if alphas[i] > 0.0:support_x.append(data_mat[i,0])support_y.append(data_mat[i,1])
for i in range(100):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(3,6)
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")
+++++++++++刚好在这一刻，2018年的蓝色月全食发生了+++++++++++++++
++++++向月全食致敬，祝我一年好运气，2018.1.31，19：41月全食++++++++
机器学习-支持向量机(python3代码实现)相关推荐
- 机器学习实战 支持向量机SVM 代码解析
机器学习实战 支持向量机SVM 代码解析 <机器学习实战>用代码实现了算法,理解源代码更有助于我们掌握算法,但是比较适合有一定基础的小伙伴.svm这章代码看起来风轻云淡,实则对于新手来说有 ...
- 机器学习 文本分类 代码_无需担心机器学习-如何在少于10行代码中对文本进行分类
机器学习 文本分类 代码 This article builds upon my previous two articles where I share some tips on how to get ...
- [系统安全] 三十三.恶意代码检测(3)基于机器学习的恶意代码检测技术
您可能之前看到过我写的类似文章,为什么还要重复撰写呢?只是想更好地帮助初学者了解病毒逆向分析和系统安全,更加成体系且不破坏之前的系列.因此,我重新开设了这个专栏,准备系统整理和深入学习系统安全.逆向分 ...
- 机器学习手撕代码（5）svm
机器学习手撕代码(5)svm 本篇分享一下svm的代码,svm.py为支持向量机模型的代码.utils.py中为可视化结果的工具. dataset见本系列第0篇. svm.py import nump ...
- 机器学习 —— 支持向量机简单入门
机器学习 -- 支持向量机简单入门 第1关:线性可分支持向量机 1.线性二分类问题 2.基本思想 3.间隔与支持向量 4.对偶问题 5.选择题 第2关:线性支持向量机 0.数据集介绍 1.软间隔 2. ...
- Python2代码转换成Python3代码
Python2代码转换成Python3代码 找到 2to3.py ,一般python安装的都是默认位置的话,位置就在: 文件夹地址栏输入: C:\Users\Administrator\AppData ...
- python3小游戏源代码_如何用python3代码玩小游戏？
在大家的印象中,程序员似乎一直是在码代码的,做着枯燥无聊的生活,殊不知,他们其实也有很多在编程中的快乐.小编最羡慕的就是他们能写一段小程序运行出来,好玩又好看,看起来还很高大上!为了照顾众多pytho ...
- 新书首发 | 《机器学习 公式推导与代码实现》正式出版！（文末送书）
大家好!我是louwill. 经过一年零三个月的努力,<机器学习 公式推导与代码实现>已于日前正式出版了. 关注过这本书的公众号读者应该知道,这本书在系列原创机器学习30讲的基础上,并参考 ...
- 温州大学《机器学习》课程代码（二）（回归）
温州大学<机器学习>课程代码(二)(回归) 代码修改并注释:黄海广,haiguang2000@wzu.edu.cn 课件 视频 下载地址:https://github.com/feng ...
最新文章
- Python的最佳学习方式
- 认识下PHP如何使用 phpmailer 发送电子邮件
- php unpack linux,PHP unpack()函数中断处理信息泄露漏洞
- Mysql数据库按照varchar字符串类型排序和按照int整型类型排序的区别和注意点及解决方案...
- telegram 机器人_学习使用Python在Telegram中构建您的第一个机器人
- 腾讯QQ团队开源分布式后台毫秒服务引擎全解析：引擎架构、RPC、灰度……
- [osg]osgDB的加载机制，使用3DS插件做参考（转，整理现有osgDB资料）
- C# List最大值最小值问题 List排序问题 List Max/Min
- matlab学习网站
- asp探针,php探针,jsp探针
- 参加2022 年第四届齐鲁工业大学（山东省科学院）与山东师范大学ICPC 大学生程序设计竞赛的总结
- 每日一题：16. “气球” 的最大数量 （C++）
- 慧安金科黄铃：面对金融欺诈， AI 如何揪出“老赖”
- 山寨手机的操作系统（mtk）简介
- 如何对磁盘分区进行写保护
- Bsp开发的几个层次
- nginx一篇入门：安装、静态网站部署、反向代理、负载均衡
- 2018年广东工业大学文远知行杯新生程序设计竞赛 1007 活在无尽梦境的后续 β
- c语言scanf返回值
- 【软件测试】我们测试人搭上元宇宙的列车，测试一直在进军......