经典SVM之SMO算法实现

                                                         经典SVM之SMO算法实现


一、浅谈理论

(1)原始问题部分

       对于理论不做过深的解释和讨论。原因有两个。第一,我也不会。第二,写公式符号太头疼!所以,只是简单给出一些最重要的公式或者程序直接使用的公式并做适当的解释。现在就给出SVM最核心最重要的一个公式:
线性支持向量机原始最优化问题:


非线性支持向量机学习算法:
【1】选取适当的核函数K(x,z)和适当的参数C,构造并求解最优化问题
           (1.1)
                          (1.2)
求得最优解
                                                                      (1.3)

【2】选择的一个正分量0 < αj < C,计算
                                        (1.4)
【3】构造决策函数:
            (1.5)
当K(x,z)是正定核函数时,(1.1)~(1.5)是凸二次规划问题,解是存在的。

(2)问题求解部分

  对于凸二次规划问题的求解有很多方法,比较流行的一种方法是用SMO(sequence minimal optimization)算法。

【1】SMO算法介绍
SMO算法是一种启发式算法,其基本思路是:如果所有变量的解都满足此最优化问题的KKT条件(1.1~1.5),那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件。否则选择两个变量,固定其他变量针对这两个变量构建一个二次规划问题。
【2】两个变量二次规划的求解方法
处理边界L,H的值:

                                        (1.6)
   时,
                                 (1.7)          
求解决策误差Ei:
                                                   (1.8)
        当i=1,2时,Ei为函数g(x)对输入xi的预测值与真实值yi之差。

对于求解未剪辑前的第二个变量的公式:
                                         (1.9)
经剪辑后的的解是:
                          (2.0)
求得
                                     (2.1)

【3】计算阀值b
               (2.2)
 
【3】总结分析

 (1)第1个变量的选择
          SMO称选择第一个变量的过程为外层循环,外层循环在训练样本中选取违反KKT条件最严重的样本点,并将其对应的变量作为第1个变量。具体地,检验训练数据样本点(xi, yi)是否满足KKT条件,即
sigmai = 0  <==>yi * g(xi) >= 1
0 < sigmai < C   <==> yi * g(xi) = 1
sigmai = C <==> yi * g(xi) <= 1
其中,g(xi) = ∑[j=1,..,N] sigmaj * yj * K(xi, xj) + b
在检验过程中,外层循环首先遍历所有满足条件0 < sigmai < C的样本点,即在间隔边界上的支持向量点,检验他们是否满足KKT条件。如果这些样本点都满足KKT条件,那么满足遍历整个训练集,检验它们是否满足KKT条件。
(2)、第2个变量的选择
SMO称选择第2个变量的过程为内层循环。假设在外层循环中已经找到了第1个变量sigma1,现在要在内层循环中找到第2个变量sigma2.第2个变量选择的标准是希望能使sigma2有最够大的变化。
sigma2[new] 是依赖|E1 - E2|的,为了加快计算速度,一种简单地做法是选择sigma2,使其对应的|E1 - E2|最大。

还有一个公式,在线性核时候会用到。计算w*
                                                                                                       (2.3)

二、代码实现
(1)数据结构
class optStruct:
	#@dataMatIn: 数据集,type: mat
	#@classLabels: 类标签,type: mat
	#@C:自设调节参数
	#@toler: 自设容错大小
	#@kTup: 核函数类型
	def __init__(self, dataMatIn, classLabels, C, toler, kTup):
		self.X = dataMatIn
		self.labelMat = classLabels
		self.C = C
		self.toler = toler
		self.m = shape(dataMatIn)[0]
		self.alphas = mat(zeros((self.m, 1)))   #初始化一个m的列向量α
		self.b = 0
		self.eCache = mat(zeros((self.m, 2)))	#误差(Ei)缓存
		self.K = mat(zeros((self.m, self.m)))   #初始化一个存储核函数值得m*m维的K
		for i in range(self.m):					#获得核函数的值K(X,xi)
			self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)


(2)核转换函数(高斯核)
#@kTup: lin: 线性 
#       rbf: 高斯核 公式:exp{(xj - xi)^2 / (2 * δ^2)} | j = 1,...,N
#       δ:有用户自设给出kTup[1]   
def kernelTrans(X, A, kTup):
	m, n = shape(X)
	K = mat(zeros(m, 1))
	if kTup[0] == 'lin': K = X * A.T
	elif kTup[0] == 'rbf': 
		for j in range(m):
			deltaRow = X[j, :] - A
			K[j] = exp(K / (-1 * kTup[1] ** 2))
	else: raise NameError('Houston we have a promble -- That kernel is not recoginzed')
	return K


(3)计算误差值Ei
#根据公式【4】,计算出Ei
def calcEk(oS, k):
	fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
	Ek = fXk - float(oS.labelMat[k])
	return Ek

(4)根据(1.9)选择第二个变量α
#随机选择第一个α
def selectJrand(i, m):
	j = i
	while (j == i):
		j = int(random.uniform(0, m))
	return j
	
#根据最大步长选择第二个变量	
def selectJ(i, oS, Ei):
	maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1, Ei]        		 #Ei保存的两维中,第一维表示Ei的有效性,第二维才是真实值
	validEcachelist = nonzero(oS.eCache[:, 0].A)[0]
	if (len(validEcachelist)) > 1:
		for k in validEcachelist:
			if k == i: continue
			Ek = calcEk(oS, k)
			deltaE = abs(Ei - Ek)
			if (deltaE > maxDeltaE):      #选择具有最大步长(max|Ei - Ej|)的j,
				maxK = k; maxDeltaE = deltaE; Ej = Ek
		return maxK, Ej
	else:                                 #随机选择第一个α
		j = selectJrand(i, oS.m)
		Ej = calcEk(oS, j)
	    return j, Ej

(5)跟新Ei值
def updatEk(oS, k):
	Ek = calcEk(oS, k)
	oS.eCache[k] = [1, Ek]

 (6)Platt SMO算法中的优化例程(选择第一个变量α)
  
#调整大于H或小于L的α值
def clipAlpha(aj, H, L):
	if aj > H:
		aj = H
	if aj < L:
	    aj = L
	return aj
	
#选择第一个变量α
def innerL(i, oS):
	Ei = calcEk(oS, i)
	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)):
		j, Ej = selectJ(i, oS, Ei)      									 #第二个α选择的启发式方法
		alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy()  	 #进行深复制处理
		
		#根据公式[6]
		if (oS.labelMat[i] != oS.labelMat[j])	:
			L = max(0, oS.alphas[j] - oS.alphas[i])
			H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
		else:
			L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
			H = min(oS.C, oS.alphas[j] + oS.alphas[i])
			
			
		if L == H: print "L == H"; return 0
		
		#根据公式[5]
		eta = 2.0 *  oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
		if eta >= 0; print "eta >= 0"; return 0                             #??为什么eta>=0就可表示剪辑过?
		oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
		oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
		updatEk(oS, j)														#跟新误差缓存
		if (abs(oS.alphas[j] - alphaJold) < 0.00001):
			print "j not moving enough"; return 0
		oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
		updatEk(oS, i)
		
		#根据公式[3]
		b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (oS.labelMat[j] - alphaJold) * oS.K[i, j]
		b2 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * (oS.labelMat[j] - alphaJold) * oS.K[j, j]
	    
		#根据SMO中优化的部分选择b的公式
		if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
		elif (o < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
		else: oS.b = (b1 + b2) / 2.0
		return 1
	else: return 0

(7)Platt SMO的外循环代码
#@maxIter:程序最大的循环次数
#其他变量的含义已经在数据结构中介绍了,这里就不在赘诉
	
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin', ):
	oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
	iter = 0
	entireSet = True; alphaPairsChanged = 0
	while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
		alphaPairsChanged = 0
		if entireSet:															#遍历所有值
			for i in range(oS.m):
				alphaPairsChanged += innerL(i, oS)
			print "fullSet, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged)
			#iter += 1
		else:																	#遍历非边界值(任意一对α值发生改变,那么就会返回1)
			nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]      #
			for i in nonBoundIs:
				alphaPairsChanged += innerL(i, oS)
				print "non-bound, iter: %d i: %d, pairs changed %d" %(iter, i, alphaPairsChanged)
			#iter += 1
		iter += 1
		if entireSet: entireSet = False
		elif (alphaPairsChanged == 0): entireSet = True
		print "iteration number: %d" % iter
	return oS.b, oS.alphas


(7)计算w值
  
#@param: alphas, dataAr, classLabels --> Type: List		
def calWs(alphas, dataArr, classLabels):
	X = mat(dataArr); labelMat = mat(classLabels).transpose()
	m, n = shape(X)
	w = zeros((n, 1))
	for i in range(m):
		w += multiply(alphas[i] * labelMat[i], X[i, :].T)
	return w

(8)根据决策函数验证结果
def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelMat[svInd];
    print "there are %d Support Vectors" % shape(sVs)[0]
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)  











本页内容版权归属为原作者,如有侵犯您的权益,请通知我们删除。
虚拟机迁移使资源配置更加灵活,尤其是在线迁移,提供了虚拟机的可用性和可靠性。Openstack liberty中提供了两种类型的迁移实现:静态迁移(cold migration)和动态迁移(live migration)。在接下来的几篇文章中,我将详细分析两种迁移的实现过程,先来看静态迁移。 限于篇幅,静态迁移的源码分析将包含两篇文章: 第一篇:主要介绍迁移过程中 nova-api 及 nova-conductor 所在的工作 第二篇:重点介绍 nova-compute 的处理过程 下面请看第一篇的内容:
  hadoop的单机模式和伪分布式模式可以说是学习hadoop的入门搭建环境,主要通过简单环境的搭建,对hadoop的MapReduce和HDFS有一个基础的认识。为分布式集群的搭建和学习起到引领的作用。   对于搭建所需的平台和软件如下:Ubuntu16.04、Hadoop2.7.2、java、sshd。以上软件都是到目前为止最新的版本。以下由于只是单机学习的目的,不考虑系统使用上的安全性,因此直接使用root用户进行操作,不再新建用户。 一、安装   1、Ubuntu16.04安装:我们通常使用虚拟

FastDFS合并存储策略 - 2016-07-25 14:07:35

FastDFS提供了合并存储功能的实现,所有的配置都在tracker.conf文件之中,具体摘录如下: trunk功能启动与配置:通过tracker.conf文件启动与配置,个配置项如下: use_trunk_file = false#是否启用trunk存储 slot_min_size = 256#trunk文件最小分配单元 slot_max_size = 16MB#trunk内部存储的最大文件,超过该值会被独立存储 trunk_file_size = 64MB#trunk文件大小 trunk_creat
论机器人的环境感知与智主运动 –兼谈基于微分几何的人工智能 标签(空格分隔): 人工智能 计算机视觉 自主移动 微分流形 Ricci流 版权声明:本文为作者原创文章,未经作者允许不得转载。 前言 人工智能是分主观与客观的,是硬币的两个方面, 客观智能是世界的本质描述,是物理的是数学的, 主观智能是来自于客观智能,是哲学的是宗教的。 抛开物理与数学只讲方法是走不远的,如无本之木、无源之水, 单讲物理与数学只会得到一个静默纷扰的世界,无乐无诗无书无画。 序言 什么是智能,这是一个令人思绪飞扬的问题,本文的内容

理解镜像、容器和存储驱动 - 2016-07-23 19:07:52

理解镜像、容器和存储驱动 为了更有效地使用存储驱动,你必须理解Docker如何创建和存储镜像。接下来,需要理解容器是如何使用这些镜像的。最后,你需要一个对镜像和容器操作者都需要的技术简介。   镜像和图层layers 每一个Docker镜像都参考了一系列的只读层,这些层代表着文件系统的区别。层级是从底层开始,逐一建立组成容器的root文件系统。下面的图显示了Ubuntu镜像有4层: Docker的存储驱动是负责堆放这些层级并且提供一个统一的视图。 当你创建一个新的容器,你会在底层栈上加入一个新的、稀疏的、
Flume简介 Flume安装 1. 解压 flume安装包到 /itcast/ 目录下 tar -zxvf /*flume安装包*/ /itcast/ 2. 修改 flume配置文件: 2.1 flume-env.sh 修改文件名称: mv flume-env.sh.template flume-env.sh 添加 java_home ,保证 flume 所使用的jdk和hdfs是一样的(可以使用 echo JAVA_HOME 查看当前机器所使用的javaHome所在路径) 2.2 编写agent配置文
注:本文为自用,随时更新。 一.系统环境 Windows7、Ecplise 4.4.0,Hadoop2.7.2 Hadoop安装教程请看: Hadoop2.7.2安装教程 此处需要注意的是:你的hadoop运行在远程的虚拟机上,但是你的windows上也需要有hadoop的运行环境,这样eclipse才能进行远程调试,所按照上面教程安装完虚拟机上的hadoop,还需要下载同样版本的Hadoop加压到windows上即可,还需要配置相应的环境变量,拷贝winutil.exe.hadoop.dll等文件到ha
2、spark wordCount程序深度剖析 标签: spark 一、Eclipse(scala IDE)开发local和cluster (一). 配置开发环境 要在本地安装好java和scala。  由于spark1.6需要scala 2.10.X版本的。推荐 2.10.4,java版本最好是1.8。所以提前我们要需要安装好java和scala并在环境变量中配置好。 下载scala IDE for eclipse安装 连接: http://scala-ide.org/download/sdk.html
1、spark 部署 标签: spark 0 apache spark项目架构 spark SQL -- spark streaming -- MLlib -- GraphX 0.1 hadoop快速搭建,主要利用hdfs存储框架 下载hadoop-2.6.0,解压,到etc/hadoop/目录下 0.2 快速配置文件 cat core-site.xml configuration property name fs.defaultFS /name value hdfs://worker1:9000 /va

Devstack单节点环境实战配置 - 2016-07-23 14:07:05

本实验是在VMware12下建立虚机的一个测试环境。 1 前期准备工作 真机环境win10 Linux版本 centos-everything-7.0 VMware版本 VMwareworkstations12 虚机配置如下: 8G内存 2核cpu(开启虚拟化) 网络配置为桥接模式 /boot 500M(一定要分大点不然之后会遇到问题) swap分区4G 其余的空间全部分给/分区 配置yum源,你可以保持装机自带的centos自带官方yum源,或者使用国内的给的镜像,本次试验中用的国外镜像并且使用fast