SVM支持向量机推导,工具介绍及python实现

支持向量机整理


参考:
Alexandre KOWALCZYK大神的SVM Tutorial
http://blog.csdn.net/alvine008/article/details/9097111
http://blog.csdn.net/zouxy09/article/details/17292011
http://blog.csdn.net/zy_zhengyang/article/details/45009431


介绍整理了SVM的基本数学推导,SMO算法的基本过程,LibSVM的用法,SMO算法的python实现四大部分(pdf版及测试数据相关电子书百度云 密码:itih)

基本型

支持向量机(Support vector machine)通常用来解决二分类问题,首先考虑一种最简单的情况,即数据线性可分的情况。

  • 划分超平面

考虑下面两点的划分,对于平面上的点集,存在一条直线将其完全划分,但是在实际情况中,样本属性通常是多维的,所以称划分线为划分超平面。

通常将划分超平面表示为:

\[\vec w·\vec x+b=0 \]

注:有的地方将\(\vec w\)表示为\(\vec w^T\) 是考虑到w和x相乘需要列数和行数相等,但是参照线性代数里可省略转置

  • Margin及支持向量

显然存在最优的划分使两类数据集距离最远,此时的泛化性能最好,称最大距离为margin,最优的划分为最优超平面,两边数据集的边界为支持向量。

首先定义支持向量为:

\[\vec w\vec x+b=s \quad or \quad \vec w\vec x+b=-s \]

同除以s可得:

\[\vec w\vec x+b=1 \quad or \quad \vec w\vec x+b=-1 \]

则对于右上方的点集y定义为1,左下点集y定义为-1

显然右上方的点集为:\(\vec w\vec x+b \ge1 \quad y=1\),左下方的点集可表示为:\(\vec w\vec x+b1\le-1 \quad y=-1\)

可合并为

\[y(\vec w\vec x+b) \ge1 \]

  • Margin的表示

计算margin如下图,长度用k表示,\(H_1\)\(H_0\)为支持向量

对于\(H_1\)上的点\(z_0\)满足

\[\vec w\vec z_0+b=1 \]

显然有

\[\vec z_0=\vec x_0 +\vec k \]

带入化简可得

\[\vec w\vec x_0+b+\vec w \vec k=1 \]

\[\vec w \left (k\frac {\vec w}{||\vec w||} \right )=2 \]

\[k=\frac {2}{||\vec w||} \]

最大的k值即为所求

  • 有约束的最优化问题

通过上面的计算,可将问题表示为约束条件下的最优化问题

\[max\quad \frac {2}{||\vec w||} \\ s.t.\quad y_i(\vec w\vec x_i+b) \ge1 \quad i=1...m \]

目标函数经常表示为 \(min\ \frac{1}{2}||\vec w||^2\)

对偶问题

通常引入拉格朗日函数简化上述问题,称此类问题为对偶问题(The wolfe dual problem)

  • 拉格朗日乘数法

拉格朗日乘数法用来解决等式约束下的最优化问题

将目标函数表示为\(f(\vec w)= \frac {2}{||\vec w||}\) ,约束函数表示为\(g(\vec w,b)=y_i(\vec w\vec x_i+b) \ge1\) ,定义如下拉格朗日函数\(L(\vec w,b,\alpha)\),其中\(\alpha_i\)为拉格朗日算子

\[L(\vec w,b,\alpha)=f(\vec w)-\sum_{i=1}^m \alpha_ig_i(\vec w,b) \\ L(\vec w,b,\alpha)=\frac{1}{2}||\vec w||^2-\sum_{i=1}^m \alpha_ig_i[y_i(\vec w\vec x_i+b)-1] \]

求解要求\(\nabla L=0\),可以解得

\[\vec w=\sum_{i=1}^m \alpha_i y_i \vec x_i \\ \sum_{i=1}^m \alpha_i y_i=0 \]

\(\vec w\)带入\(L(\vec w,b,\alpha)\)可得

\[f(\vec w,b)=\sum_{i=1}^m\alpha_i-\frac{1}{2} \sum_{i=1}^m\sum_{j=1}^m \alpha_i\alpha_jy_iy_j \vec x_i \vec x_j \]

所以问题转化为

\[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 \vec x_i \vec x_j \\ s.t. \quad \alpha_i\ge0 \quad i=1...m \\ \quad \sum_{i=1}^m \alpha_i y_i=0 \]

  • KKT条件

注意到上面的拉格朗日乘数法用来解决等式优化问题,在对于不等式约束的优化问题时还需附加如下的KKT条件

\[\begin{cases} \ \alpha_i \ge 0 \\ \ y_if(x_i)-1 \ge 0 \\ \ \alpha_i(y_if(\vec x_i)-1) \ge 0 \end{cases} \]

软间隔

  • 线性不可分时

实际遇到的问题数据总是有噪音的,也就是说可能存在一两个可忽略的特殊点使margin大大减小,也可能有一两个特殊点导致集合线性不可分(如下图),所以我们需要允许模型犯一定的错误来过滤掉噪音,即允许间隔的调整称为软间隔

  • 松弛变量

引入松弛变量(slack variable)\(\xi_i\) ,则约束函数可写为

\[y_i(\vec w\vec x_i+b) \ge1-\xi_i \]

  • 控制出错

同时引入变量C表示对\(\xi_i\)的控制程度,即表示允许模型忽略噪音的程度,显然当C->\(\infty\)退化成基本型的硬间隔,C->0时成为无约束问题

引入C后基本型可表示为

\[max\quad \frac {1}{2}{||\vec w||}^2+C\sum_{i=1}^m \xi_i \\ s.t.\quad y_i(\vec w\vec x_i+b) \ge1-\xi_i \quad \\ \xi_i \ge 0 \quad i=1...m \]

同样由拉格朗日乘数法可得

\[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\vec x_i \vec x_j \\ s.t. \quad 0\le \alpha_i\le C \quad i=1...m \\ \quad \sum_{i=1}^m \alpha_i y_i=0 \]

核方法

  • 升维

考虑更通常的情况,当点集线性不可分时,更有趣的方法是升维,如下图在二维平面内不能找到一个划分超平面

然而当在更高维的情况比如三维,我们就可以找到划分超平面来解决问题

这种方法通常成为核方法(kernel)

  • 构造核函数

比如我们可以简单的\((\vec x,\vec y)\) ->\((\vec x,\vec x\vec y,\vec y)\) 将点集升维,但是当样本空间很大时代价太高,考虑上面的目标函数

我们可以利用目标函数中的点积来构建核函数避免计算量的增加

常用核函数

  1. 线性核:\(k(\vec x_i,\vec y_i)=x_i·y_i\)
  2. 多项式核:\(k(\vec x,\vec y)=(\vec x_i·\vec y_i)^d\)\(d\ge 1\)为多项式的次数
  3. 高斯核(RBF,常用):\(k(\vec x,\vec y)=exp(-\frac {||x_i-x_j||^2}{2 \sigma ^2})\)\(\sigma\)为高斯核的带宽

SMO算法

  • SMO概念

对于上述问题的求解,SMO算法每次迭代优化两个\(\alpha\) (若只选择1个,由约束条件可以直接解出\(\alpha\),就不是最优化问题了),把原始求解N个参数二次规划问题分解成很多个子二次规划问题分别求解,每个子问题只需要求解2个参数,方法类似于坐标上升,节省时间成本和降低了内存需求。每次启发式选择两个变量进行优化,不断循环,直到达到函数最优值。

  1. 视为一个二元函数

SMO算法选择同时优化两个参数,固定其他N-2个参数,假设选择的变量为α1,α2,固定其他参数α3,α4,...,αN,由于参数α3,α4,...,αN的固定,可以简化目标函数为只关于α1,α2的二元函数,Const表示常数项(不包含变量α1,α2的项)。

\[minW(\alpha_i,\alpha_j)=\frac {1}{2}K_{11}\alpha_1^2+\frac {1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1 \alpha_2-(\alpha_1+\alpha_2)+y_1v_1\alpha_1+y_2v_2\alpha_2+Const \\ 其中,v_i=\sum _{j=3}^m \alpha_jy_jK(x_i,x_j) \quad i=1,2 \]

  1. 视为一元函数求极值点

由等式约束:\(\alpha_i y1+\alpha_2 y2=-\sum _{i=3}^{m} \alpha_i y_i=\xi\),可见\(\xi\)为定值,将\(\alpha_1 y1+\alpha_2 y2=\xi\)两边同时乘以\(y_1\),且\(y_1^2=1\),得

\[\alpha_1=(\xi-y_2\alpha_2)y_1 \]

代回上述二元函数可得关于\(\alpha_2\)的一元函数,求导令其等于0可得

\[\frac {\partial W(\alpha_2)}{\partial \alpha_2}=(K_{11}+K_{22}-2K_{12})\alpha_2-K_{11}\xi y_2+K_{12}\xi y_2+y_1y_2-1-v_1y_1+v_2y_2=0 \]

假设从上式中求得了\(\alpha_2\)的解,代回可得\(\alpha_1\)分别记为\(\alpha_1^{new} ,\alpha_2^{new}\)优化前的为\(\alpha_1^{old} ,\alpha_2^{old}\)

同样由约束条件可得

\[\alpha_1^{old} y1+\alpha_2^{old} y2=\xi \]

在对偶问题中解得的\(\vec w\)代入超平面方程\(f(x)=\vec w\vec x+b=\sum_{i=1}^{m}\alpha_iy_iK(x_i,x_j)+b\),同时定义\(E_i\)表示预测值与真实值之差

\[E_i=f(x_i)-y_i \]

\(v_i=\sum _{j=3}^m \alpha_jy_jK(x_i,x_j) \quad i=1,2\)可得

\[v_1=f(x_1)-\sum_{j=1}^{2}y_j\alpha_jK_{1j}-b \\ v_2=f(x_2)-\sum_{j=1}^{2}y_j\alpha_jK_{2j}-b \]

将上述式子带入\(\frac {\partial W(\alpha_2)}{\partial \alpha_2}=0\)可得

\[\alpha_2^{new}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta},其中\eta=K_{11}+K{22}-2K_{12} \]

  1. 对上述解进行修剪

考虑到约束条件

\[0\le \alpha_i \le C \\ \alpha_1y_1+\alpha_2y_2=\xi \]

定义\(\alpha_2\)区间为[L,H],则有

\(y_1\not=y_2\)时,\(L=max(0,\alpha_2^{old}-\alpha_1^{old})\);\(H=min(C,C+\alpha_2^{old}-\alpha_1^{old})\);

\(y_1=y_2\)时,\(L=max(0,\alpha_2^{old}+\alpha_1^{old}-C)\);\(H=min(C,\alpha_2^{old}+\alpha_1^{old})\);

  1. 求解\(\alpha_1^{new}\)

\(\alpha_1^{old} y1+\alpha_2^{old}y2=\alpha_1^{new} y1+\alpha_2^{new}y2\) ,可得

\[\alpha_1^{new}=\alpha_1^{old} +y_1y_2(\alpha_2^{old}-\alpha_2^{new}) \]

  1. 考虑临界情况

大部分情况下由\(\eta=K_{11}+K_{22}-2K_{12}>0\),但当如下两种情况时\(\alpha_2\)取临界值L或H

  • η<0,当核函数K不满足Mercer定理时,矩阵K非正定;
  • η=0,样本x1与x2输入特征相同;

也可以理解成对一元函数求二阶导就是\(\eta=K_{11}+K_{22}-2K_{12}>0\)

  • \(\eta\)<0时,目标函数为凸函数,没有极小值,极值在定义域边界处取得
  • \(\eta\)=0时,目标函数为单调函数,同样在边界处取极值
  1. 启发式选择变量

上述分析是在从N个变量中已经选出两个变量进行优化的方法,下面分析如何高效地选择两个变量进行优化,使得目标函数下降的最快。

  • \(\alpha_1\)

第一个变量的选择称为外循环

先遍历非边界样本集(0<\(\alpha_1\)<C)中违反KKT的\(\alpha_1\)作为第一个变量,同样依据相关规则选择第二个变量,对此两个变量进行优化。

当遍历完非边界样本集后,再次回到边界样本集中寻找,即在边界样本集与非边界样本集上来回切换,寻找违反KKT条件的\(\alpha_1\)作为第一个变量。直到遍历整个样本集后,没有违反KKT条件αi,然后退出。
边界上的样本对应的\(\alpha_1\)=0或者\(\alpha_1\)=C,在优化过程中很难变化,然而非边界样本0<\(\alpha_1\)<C会随着对其他变量的优化会有大的变化。

  • \(\alpha_2\)

第二个变量的选择过程为内循环

第二个变量的选择希望能使\(\alpha_2\)有较大的变化,通常为每个样本的\(E_i\)保存在一个列表中,选择最大的\(|E_1-E_2|\)来近似最大化步长,由于α2是依赖于\(|E_1-E_2|\),当\(E_1\)为正时,那么选择最小的\(E_i\)作为\(E_2\),如果\(E_1\)为负,选择最大\(E_i\)作为\(E_2\)

  1. 更新b

每完成对两个变量的优化后,要对b的值进行更新,因为b的值关系到f(x)的计算,即关系到下次优化时\(E_I\)的计算。

LibSVM工具包求解

LibSVM是台湾林智仁(Chih-Jen Lin)教授2001年开发的一套支持向量机的库,这套库运算速度快,可以很方便的对数据做分类或回归。由于libSVM程序小,运用灵活,输入参数少,并且是开源的,易于扩展,因此成为目前国内应用最多的SVM的库。(以下参考博客)

  • 基本介绍

可以从https://www.csie.ntu.edu.tw/~cjlin/libsvm/下载,解压到任意目录,有如下几个目录

Java——主要是应用于java平台;

Python——是用来参数优选的工具,稍后介绍;

svm-toy——一个可视化的工具,用来展示训练数据和分类界面,里面是源码,其编译后的程序在windows文件夹下;

tools——主要包含四个python文件,用来数据集抽样(subset),参数优选(grid),集成测试(easy),数据检查(checkdata);

windows——包含libSVM四个exe程序包,一般直接用这四个就够了

其他.h和.cpp文件都是程序的源码,可以编译出相应的.exe文件。

  • 自带数据测试
  1. 命令行下进入libsvm下的windows目录
  2. 输入命令 .\svm-train ..\heart_scale test_model

heart_scale是目录下的已经存在的样本文件,以后改成自己的文件名就可以了, test_model——是创建模型

其中,#iter为迭代次数,

nu是你选择的核函数类型的参数

obj为SVM文件转换为的二次规划求解得到的最小值

rho为判决函数的偏置项b

nSV为标准支持向量个数(0<a[i]<c)

nBSV为边界上的支持向量个数(a[i]=c),

Total nSV为支持向量总个数(对于二分类来说,因为只有一个分类模型Total nSV = nSV,但是对于多类,这个是各个分类模型的nSV之和)。

在目录下,还可以看到产生了一个train.model文件,可以用记事本打开,记录了训练后的结果。

svm_type c_svc//所选择的svm类型,默认为c_svc

kernel_type rbf//训练采用的核函数类型,此处为RBF核

gamma 0.0769231//RBF核的参数γ

nr_class 2//类别数,此处为两分类问题

total_sv 132//支持向量总个数

rho 0.424462//判决函数的偏置项b

label 1 -1//原始文件中的类别标识

nr_sv 64 68//每个类的支持向量机的个数

SV//以下为各个类的权系数及相应的支持向量

  • libsvm使用规范
  1. libsvm的数据格式

Label 1:value 2:value(value为0时序号可直接省略)

  1. svm-train的用法

svmtrain我们在前面已经接触过,他主要实现对训练数据集的训练,并可以获得SVM模型。

用法: svmtrain [options] training_set_file [model_file]

其中,options为操作参数,可用的选项即表示的涵义如下所示:

-s设置svm类型:

0 – C-SVC

1 – v-SVC

2 – one-class-SVM

3 –ε-SVR

4 – n - SVR

-t设置核函数类型,默认值为2

0 --线性核:u'**v*

1 --多项式核:(gu'*v+coef0)degree*

2 -- RBF核:exp(-γ*||u-v||2)

3 -- sigmoid核:tanh(γu'*v+coef*0)

-d degree:设置多项式核中degree的值,默认为3

-gγ:设置核函数中γ的值,默认为1/k,k为特征(或者说是属性)数;

-r coef 0:设置核函数中的coef 0,默认值为0;

-c cost:设置C-SVC、ε-SVR、n - SVR中从惩罚系数C,默认值为1;

-n v:设置v-SVC、one-class-SVM与n - SVR中参数n,默认值0.5;

-pε:设置v-SVR的损失函数中的e,默认值为0.1;

-m cachesize:设置cache内存大小,以MB为单位,默认值为40;

-eε:设置终止准则中的可容忍偏差,默认值为0.001;

-h shrinking:是否使用启发式,可选值为0或1,默认值为1;

-b概率估计:是否计算SVC或SVR的概率估计,可选值0或1,默认0;

-wi weight:对各类样本的惩罚系数C加权,默认值为1;

-v n:n折交叉验证模式;

model_file:可选项,为要保存的结果文件,称为模型文件,以便在预测时使用。

默认情况下,只需要给函数提供一个样本文件名就可以了,但为了能保存结果,还是要提供一个结果文件名,比如:test.model,则命令为:svmtrain test.txt test.model

  • svm-predict的用法

svmpredict是根据训练获得的模型,对数据集合进行预测。

用法:svmpredict [options] test_file model_file output_file

其中,options为操作参数,可用的选项即表示的涵义如下所示:

-b probability_estimates——是否需要进行概率估计预测,可选值为0或者1,默认值为0。

model_file ——是由svmtrain产生的模型文件;

test_file——是要进行预测的数据文件,格式也要符合libsvm格式,即使不知道label的值,也要任意填一个,svmpredict会在output_file中给出正确的label结果,如果知道label的值,就会输出正确率;

output_file ——是svmpredict的输出文件,表示预测的结果值。

Python实现SMO算法

#smo
from numpy import *
import time


# calulate kernel value
def calcKernelValue(matrix_x, sample_x, kernelOption):
	kernelType = kernelOption[0]
	numSamples = matrix_x.shape[0]
	kernelValue = mat(zeros((numSamples, 1)))

	if kernelType == 'linear':
		kernelValue = matrix_x * sample_x.T
	elif kernelType == 'rbf':
		sigma = kernelOption[1]
		if sigma == 0:
			sigma = 1.0
		for i in xrange(numSamples):
			diff = matrix_x[i, :] - sample_x
			kernelValue[i] = exp(diff * diff.T / (-2.0 * sigma**2))
	else:
		raise NameError('Not support kernel type! You can use linear or rbf!')
	return kernelValue


# calculate kernel matrix given train set and kernel type
def calcKernelMatrix(train_x, kernelOption):
	numSamples = train_x.shape[0]
	kernelMatrix = mat(zeros((numSamples, numSamples)))
	for i in xrange(numSamples):
		kernelMatrix[:, i] = calcKernelValue(train_x, train_x[i, :], kernelOption)
	return kernelMatrix


# define a struct just for storing variables and data
class SVMStruct:
	def __init__(self, dataSet, labels, C, toler, kernelOption):
		self.train_x = dataSet # each row stands for a sample
		self.train_y = labels  # corresponding label
		self.C = C             # slack variable
		self.toler = toler     # termination condition for iteration
		self.numSamples = dataSet.shape[0] # number of samples
		self.alphas = mat(zeros((self.numSamples, 1))) # Lagrange factors for all samples
		self.b = 0
		self.errorCache = mat(zeros((self.numSamples, 2)))
		self.kernelOpt = kernelOption
		self.kernelMat = calcKernelMatrix(self.train_x, self.kernelOpt)


# calculate the error for alpha k
def calcError(svm, alpha_k):
	output_k = float(multiply(svm.alphas, svm.train_y).T * svm.kernelMat[:, alpha_k] + svm.b)
	error_k = output_k - float(svm.train_y[alpha_k])
	return error_k


# update the error cache for alpha k after optimize alpha k
def updateError(svm, alpha_k):
	error = calcError(svm, alpha_k)
	svm.errorCache[alpha_k] = [1, error]


# select alpha j which has the biggest step
def selectAlpha_j(svm, alpha_i, error_i):
	svm.errorCache[alpha_i] = [1, error_i] # mark as valid(has been optimized)
	candidateAlphaList = nonzero(svm.errorCache[:, 0].A)[0] # mat.A return array
	maxStep = 0; alpha_j = 0; error_j = 0

	# find the alpha with max iterative step
	if len(candidateAlphaList) > 1:
		for alpha_k in candidateAlphaList:
			if alpha_k == alpha_i:
				continue
			error_k = calcError(svm, alpha_k)
			if abs(error_k - error_i) > maxStep:
				maxStep = abs(error_k - error_i)
				alpha_j = alpha_k
				error_j = error_k
	# if came in this loop first time, we select alpha j randomly
	else:
		alpha_j = alpha_i
		while alpha_j == alpha_i:
			alpha_j = int(random.uniform(0, svm.numSamples))
		error_j = calcError(svm, alpha_j)

	return alpha_j, error_j


# the inner loop for optimizing alpha i and alpha j
def innerLoop(svm, alpha_i):
	error_i = calcError(svm, alpha_i)

	### check and pick up the alpha who violates the KKT condition
	## satisfy KKT condition
	# 1) yi*f(i) >= 1 and alpha == 0 (outside the boundary)
	# 2) yi*f(i) == 1 and 0<alpha< C (on the boundary)
	# 3) yi*f(i) <= 1 and alpha == C (between the boundary)
	## violate KKT condition
	# because y[i]*E_i = y[i]*f(i) - y[i]^2 = y[i]*f(i) - 1, so
	# 1) if y[i]*E_i < 0, so yi*f(i) < 1, if alpha < C, violate!(alpha = C will be correct)
	# 2) if y[i]*E_i > 0, so yi*f(i) > 1, if alpha > 0, violate!(alpha = 0 will be correct)
	# 3) if y[i]*E_i = 0, so yi*f(i) = 1, it is on the boundary, needless optimized
	if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\
		(svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):

		# step 1: select alpha j
		alpha_j, error_j = selectAlpha_j(svm, alpha_i, error_i)
		alpha_i_old = svm.alphas[alpha_i].copy()
		alpha_j_old = svm.alphas[alpha_j].copy()

		# step 2: calculate the boundary L and H for alpha j
		if svm.train_y[alpha_i] != svm.train_y[alpha_j]:
			L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])
			H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])
		else:
			L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)
			H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])
		if L == H:
			return 0

		# step 3: calculate eta (the similarity of sample i and j)
		eta = 2.0 * svm.kernelMat[alpha_i, alpha_j] - svm.kernelMat[alpha_i, alpha_i] \
				  - svm.kernelMat[alpha_j, alpha_j]
		if eta >= 0:
			return 0

		# step 4: update alpha j
		svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta

		# step 5: clip alpha j
		if svm.alphas[alpha_j] > H:
			svm.alphas[alpha_j] = H
		if svm.alphas[alpha_j] < L:
			svm.alphas[alpha_j] = L

		# step 6: if alpha j not moving enough, just return
		if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:
			updateError(svm, alpha_j)
			return 0

		# step 7: update alpha i after optimizing aipha j
		svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \
								* (alpha_j_old - svm.alphas[alpha_j])

		# step 8: update threshold b
		b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
													* svm.kernelMat[alpha_i, alpha_i] \
							 - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
													* svm.kernelMat[alpha_i, alpha_j]
		b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
													* svm.kernelMat[alpha_i, alpha_j] \
							 - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
													* svm.kernelMat[alpha_j, alpha_j]
		if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):
			svm.b = b1
		elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):
			svm.b = b2
		else:
			svm.b = (b1 + b2) / 2.0

		# step 9: update error cache for alpha i, j after optimize alpha i, j and b
		updateError(svm, alpha_j)
		updateError(svm, alpha_i)

		return 1
	else:
		return 0


# the main training procedure
def trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('rbf', 1.0)):
	# calculate training time
	startTime = time.time()

	# init data struct for svm
	svm = SVMStruct(mat(train_x), mat(train_y), C, toler, kernelOption)

	# start training
	entireSet = True
	alphaPairsChanged = 0
	iterCount = 0
	# Iteration termination condition:
	# 	Condition 1: reach max iteration
	# 	Condition 2: no alpha changed after going through all samples,
	# 				 in other words, all alpha (samples) fit KKT condition
	while (iterCount < maxIter) and ((alphaPairsChanged > 0) or entireSet):
		alphaPairsChanged = 0

		# update alphas over all training examples
		if entireSet:
			for i in xrange(svm.numSamples):
				alphaPairsChanged += innerLoop(svm, i)
			print '---iter:%d entire set, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
			iterCount += 1
		# update alphas over examples where alpha is not 0 & not C (not on boundary)
		else:
			nonBoundAlphasList = nonzero((svm.alphas.A > 0) * (svm.alphas.A < svm.C))[0]
			for i in nonBoundAlphasList:
				alphaPairsChanged += innerLoop(svm, i)
			print '---iter:%d non boundary, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
			iterCount += 1

		# alternate loop over all examples and non-boundary examples
		if entireSet:
			entireSet = False
		elif alphaPairsChanged == 0:
			entireSet = True

	print 'Congratulations, training complete! Took %fs!' % (time.time() - startTime)
	return svm


# testing your trained svm model given test set
def testSVM(svm, test_x, test_y):
	test_x = mat(test_x)
	test_y = mat(test_y)
	numTestSamples = test_x.shape[0]
	supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]
	supportVectors 		= svm.train_x[supportVectorsIndex]
	supportVectorLabels = svm.train_y[supportVectorsIndex]
	supportVectorAlphas = svm.alphas[supportVectorsIndex]
	matchCount = 0
	for i in xrange(numTestSamples):
		kernelValue = calcKernelValue(supportVectors, test_x[i, :], svm.kernelOpt)
		predict = kernelValue.T * multiply(supportVectorLabels, supportVectorAlphas) + svm.b
		if sign(predict) == sign(test_y[i]):
			matchCount += 1
	accuracy = float(matchCount) / numTestSamples
	return accuracy

数据读写及输出的控制程序如下

from numpy import *
import SVM
import csv


################## test svm #####################
## step 1: load data
print "step 1: load data..."
dataSet = []
labels = []

in_file=csv.reader(open('spammer_order.csv','r'))

maxm=[0]*35
minm=[0]*35
index=1
for in_line in in_file:
    if index<3:
        if index==1:
            for i in range(4,32):
                maxm[i]=float(in_line[i])
        if index==2:
            for i in range(4,32):
                minm[i]=float(in_line[i])
    else:
        out_line=[]
        if in_line[32]=='yes':
            labels.append('1')
        else:
            labels.append('-1')

        for i in range(4,32):
            if i!=6:
                if i==4:
                    if in_line[4]=='female':
                        out_line.append('1')
                    else:
                        out_line.append('-1')
                else:
                    max_value=maxm[i]
                    min_value=minm[i]
                    value=float(in_line[i])
                    n=2*(value-min_value)/(max_value-min_value)-1
                    out_line.append(str(n))
        dataSet.append(map(float,out_line))
    index=index+1
#3 860 1718
p1=3+(860-3)/3*2
p=860
p2=860+(1718-860)/3*2
labels=map(float,labels)
train_data=mat(dataSet[0:p1]+dataSet[p:p2])
train_labels=mat(labels[0:p1]+labels[p:p2]).T
test_data=mat(dataSet[p1:p]+dataSet[p2:])
test_labels=mat(labels[p1:p]+labels[p:p2]).T

train_x = train_data
train_y = train_labels
test_x = test_data
test_y = test_labels

## step 2: training...
print "step 2: training..."
C = 1
toler = 0.001
maxIter = 428
svmClassifier = SVM.trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('rbf', 0.440787))

## step 3: testing
print "step 3: testing..."
accuracy = SVM.testSVM(svmClassifier, test_x, test_y)

## step 4: show the result
print "step 4: show the result..."
print 'The classify accuracy is: %.3f%%' % (accuracy * 100)
posted @ 2017-11-17 16:39  Neptune15  阅读(2232)  评论(0编辑  收藏  举报