SVM(支持向量机原理+代码)

最新推荐文章于 2024-09-09 16:50:48 发布

qq_35575840

最新推荐文章于 2024-09-09 16:50:48 发布

阅读量1.2k

点赞数

本文链接：https://blog.csdn.net/qq_35575840/article/details/107700534

版权

SVM

1. 支持向量机原理（硬间隔支持向量机）
2. 软间隔支持向量机
3. 非线性支持向量机（核函数）
4.多分类支持向量机
- 4.1一对多法（one-versus-rest,简称OVR SVMs）
- 4.2 一对一法（one-versus-one,简称OVO SVMs或者pairwise）
5. 支持向量回归
6. 代码
- 6.1数据
- 6.2运行

SVM，全称是support vector machine，中文名叫支持向量机。它的目标是为确定一个分类超平面，从而将不同的数据分隔开。
支持向量机学习方法包括构建由简至繁的模型,可以分为线性可分支持向量机（硬间隔最大化）、线性支持向量机（软间隔最大化）、非线性支持向量机（核技巧+软间隔最大化）。

1. 支持向量机原理（硬间隔支持向量机）

1.1 支持向量机基本型

在这里插入图片描述
在样本空间中，超平面 $(w, b)$ 可用如下方程来描述：
$\omega^Tx + b = 0 \tag{1}$
样本空间中，任意点 $x$ 到超平面到 $(w, b)$ 到距离为：
$\frac{|\omega^Tx +ｂ|}{|| \omega||} \tag{2}$
其中， $||\omega||$ 为 $\omega$ 的范数。
假设超平面 $(\omega,b)$ 能将训练样本正确分类，及对于 $(x_i,y_i) \in D$ ,若 $y_i=+1$ ,则有 $\omega^Tx_i+b>0$ ;若 $y_i=-1$ ,则有 $\omega^Tx_i+b<0$ ;令：
$\begin{Bmatrix} \omega x_i + b \geqslant +1,y_i=+1\\ \omega x_i + b \leqslant -1,y_i=-1 \tag{3} \end{Bmatrix}$
则使得公式点成立的样本点被称为“支持向量”，则两个异类支持向量到超平面的距离为：
$\gamma = \frac{2}{||\omega||} \tag{4}$
欲找到具有最大间隔的划分超平面，也就是能找到满足约束条件的参数 $\omega$ 和 $b$ ，使得 $\gamma$ 最大，即：
$\begin{matrix} \underset{\omega,b}{max} \frac{2}{||\omega||}\\ s.t. \quad y_i(\omega^Tx_i+b) \geqslant1,i=1,2,\cdots ,m. \end{matrix} \tag{5}$
所以支持向量机（SVM）的基本型可以写为：
$\begin{matrix} \underset{\omega,b}{min} \ \frac{1}{2}||\omega||\\ s.t. \quad y_i(\omega^Tx_i+b) \geqslant1, \ i=1,2,\cdots ,m. \end{matrix} \tag{6}$
因为现在的目标函数是二次的，约束条件是线性的，所以它是一个凸二次规划问题。
可以通过QP (Quadratic Programming) 优化包进行求解，也可以通过拉格朗日对偶性（Lagrange Duality）变换到对偶变量 (dual variable) 的优化问题。

1.2 对偶问题

对基本型（6）的每条约束添加拉格朗日乘子 $\alpha_i\geqslant 0$ ,则该问题的拉格朗日函数可以写为：
$L(\omega,b,\alpha) = \frac{1}{2}||\omega||^2+\sum_{i=1}^{m}\alpha_i(1-y_i(\omega^Tx_i+b)), \tag{7}$
令 $L(\omega,b,\alpha)$ 对 $\omega,b$ 偏导为零可得：
$\begin{matrix} \omega = \sum_{i=1}^m\alpha_iy_ix_i \\ \\ 0 =\sum_{i=1}^m\alpha_iy_i \end{matrix}\tag{8}$ 带入得到对偶问题：
$\begin{matrix} \underset{\alpha}{max} \quad \sum_i^m\alpha_i-\frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_jx_i^Tx_j \\ \\ s.t. \quad \sum_{i=1}^m \alpha_iy_i=0,\\ \alpha_i \geqslant 0, \quad i = 1,2,\cdots,m. \end{matrix} \tag{9}$
解出 $\alpha$ 后，求出 $\omega$ 和 $b$ 即可得到模型：
$\omega^Tx+b=\sum_{i=1}^m\alpha_iy_ix_i^Tx+b \tag{10}$

1.3 SMO算法

求解公式（9）为二次规划问题，采用SMO算法进行计算。
SMO算法的基本思想是先固定 $\alpha_i$ 之外的所有参数，然后求 $\alpha_i$ 上的极值。由于存在约束 $\sum_{i=1}^m\alpha_iy_i=0$ ,若固定 $\alpha_i$ 之外的其他变量，则 $\alpha_i$ 则可有其他变量导出。于是，SMO每次选择两个变量 $\alpha_i$ 和 $\alpha_j$ ,并固定其他参数，不断迭代更新 $\alpha_i$ 和 $\alpha_j$ 。

2. 软间隔支持向量机

在这里插入图片描述
硬间隔支持向量机要求所有样本必须满足约束（3），软间隔支持向量机允许某些样本不满足约束。
线性不可分意味着某些样本点 $x_i,y_i)$ 不能满足函数间隔大于等于1的约束条件，为了解决这个问题，可以对每个样本点 $x_i,y_i)$ 引进一个松弛变量 $\zeta_i \geqslant 0$ ，这样，约束条件变为：
$y_i( \omega x_i+b)\geqslant1-\zeta_i\tag{11}$
同时，对于每个松弛变量 $\zeta_i$ ,支付一个代价 $\zeta_i$ ,目标函数由原来的 $\frac{1}{2}||\omega||^2$ 变为：
$\frac{1}{2}||\omega||^2+C\sum_{I=1}^N\zeta_i \tag{12}$
这里，C>0称为惩罚参数，一般由应用问题决定，C值大时对误分类的惩罚增大，C值小时对误分类的惩罚减小，此时，最小化目标函数有两层含义：使 $\frac{1}{2}||\omega||^2$ 尽量小，同时使误分类的个数尽量少，C是调和二者的系数。
松弛变量 $\zeta_i$ 应是 $y_i(\omega^Tx_i+b)-1$ 的函数，即优化目标应该是：
$\underset{\omega,b}{min}\frac{1}{2}||\omega||^2+C\sum_{I=1}^N\ell_{0/1}(y_i(\omega^Tx_i+b)-1) \tag{12}$
其中， $\ell_{0/1}$ 非凸、非连续性，因此常用其他一些函数来替代：
hinge损失： $\ell_{hinge}(z)=max(0,1-z)$
指数损失： $\ell_{exp}(z)=exp(-z)$
对率损失： $\ell_{log}=log(1+exp(-z))$
之后的处理和硬间隔支持向量机的处理一样。

3. 非线性支持向量机（核函数）

硬间隔支持向量机和软间隔支持向量机都是假设样本是线性可分的或者近似线性可分的，但是现实中有很多都是线性不可分的。对于这样的问题，可以将样本从原始空间映射到一个更高维的特征空间，使得样本在这个特征空间内线性可分。
令 $\phi(x)$ 表示将 $x$ 映射后的特征向量，于是在线性空间中划分超平面所对应的模型可表示为：
$f(x)=\omega^T\phi(x)+b, \tag{13}$
其中 $\omega$ 和 $b$ 是模型参数，对比公式6有：
$\begin{matrix} \underset{\omega,b}{min} \quad \frac{1}{2}||\omega||^2\\ s.t. \quad y_i(\omega^T\phi(x_i)+b)\geqslant1, \quad i=1,2,\cdots,m. \end{matrix} \tag{14}$
其对偶问题是：
$\begin{matrix} \underset{\alpha}{max} \quad \sum_{i=1}^m\alpha_i-\frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_j\phi(x_i)^T\phi(x_j),\\ s.t. \quad \sum_{i=1}^m\alpha_iy_i=0,\\ \qquad \qquad \qquad \quad \alpha_i\geqslant0, \quad i=1,2,\cdots,m. \end{matrix} \tag{15}$
求解中，涉及 $\phi(x_i)^T\phi(x_j)$ ,这是样本 $x_i$ , $x_j$ 映射到特征空间的之后的内积，直接 $\phi(x_i)^T\phi(x_j)$ 是困难的。为了避开这个障碍，可以设想这样一个函数：
$\kappa (x_i,x_j)=<\phi(x_i),\phi(x_j)>= \phi(x_i)^T\phi(x_j) \tag{16}$
即 $x_i$ 与 $x_j$ 在特征空间的内积等于他们在原始样本空间中通过函数 $\kappa(\cdot,\cdot)$ 计算的结果。
公式（15）可以改写为：
$\begin{matrix} \underset{\alpha}{max} \quad \sum_{i=1}^m\alpha_i-\frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_j\kappa (x_i,x_j),\\ s.t. \quad \sum_{i=1}^m\alpha_iy_i=0,\\ \qquad \qquad \qquad \quad \alpha_i\geqslant0, \quad i=1,2,\cdots,m. \end{matrix} \tag{17}$
求解后可以得到：
$\begin{matrix} f(x) = \omega^T\phi(x)+b \\ \qquad \qquad \qquad \qquad =\sum_{i=1}^m\alpha_iy_i\phi(x_i)^T\phi(x)+b \\ \qquad \qquad \qquad =\sum_{i=1}^m\alpha_iy_i\kappa(x,x_i)+b \end{matrix} \tag{18}$
这里的函数 $\kappa(\cdot,\cdot)$ 就是“核函数”，公式（18）展示了模型的最优解可通过训练样本的核函数展开。
根据核函数定理:任何一个对称函数所对应的核矩阵半正定，他就能作为核函数使用。也就是，任何一个核函数都隐式的定义了一个“再生核希尔伯特空间”，几种常见的核函数如下：

名称	表达式	参数
线性核	$\kappa(x_i,x_j)=x_i^Tx_j$
多项式核	$\kappa(x_i,x_j)=（x_i^Tx_j)^d$	$\geqslant1$ 为多项式的次数
高斯核	$\kappa(x_i,x_j)=exp(- \frac{ \\|x_i-x_j\\|^2 }{2\sigma^2})$	$\sigma >0$ 为高斯核的带宽
拉普拉斯核	$\kappa(x_i,x_j)=exp(- \frac{ \\|x_i-x_j\\| }{\sigma})$	$\sigma >0$
Sigmoid核	$\kappa(x_i,x_j)=tanh(\beta x_i^Tx_j+\theta)$	tanh为双曲正切函数， $\beta>0,\theta<0$

4.多分类支持向量机

目前，构造SVM多类分类器的方法主要有两类
　　（1）直接法，直接在目标函数上进行修改，将多个分类面的参数求解合并到一个最优化问题中，通过求解该最优化问题“一次性”实现多类分类。这种方法看似简单，但其计算复杂度比较高，实现起来比较困难，只适合用于小型问题中；
　　（2）间接法，主要是通过组合多个二分类器来实现多分类器的构造，常见的方法有one-against-one和one-against-all两种。

4.1一对多法（one-versus-rest,简称OVR SVMs）

训练时依次把某个类别的样本归为一类,其他剩余的样本归为另一类，这样k个类别的样本就构造出了k个SVM。分类时将未知样本分类为具有最大分类函数值的那类。

4.2 一对一法（one-versus-one,简称OVO SVMs或者pairwise）

其做法是在任意两类样本之间设计一个SVM，因此k个类别的样本就需要设计k(k-1)/2个SVM。
当对一个未知样本进行分类时，最后得票最多的类别即为该未知样本的类别。

5. 支持向量回归

支持向量回归假设我们能容忍 $f (x)$ 与 $y$ 之间最多有 $\varepsilon$ 的偏差，即仅当 $f (x)$ 与 $y$ 之间的差别绝对值小于 $\varepsilon$ ，我们认为预测正确。
在这里插入图片描述
SVM是要使到超平面最近的样本点的“距离”最大；
SVR则是要使到超平面最远的样本点的“距离”最小。

6. 代码

6.1数据

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

6.2运行

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