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个参数，方法类似于坐标上升，节省时间成本和降低了内存需求。每次启发式选择两个变量进行优化，不断循环，直到达到函数最优值。

• SMO原理分析(参考博客)

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)