机器学习--支持向量机SMO算法1

最新推荐文章于 2023-03-23 14:08:48 发布

暗夜猎手-大魔王

最新推荐文章于 2023-03-23 14:08:48 发布

阅读量674

点赞数 1

分类专栏：机器学习机器学习

本文链接：https://blog.csdn.net/u014106644/article/details/83444409

版权

机器学习同时被 2 个专栏收录

35 篇文章 3 订阅

订阅专栏

机器学习

34 篇文章 3 订阅

订阅专栏

结合论文《Sequential Minimal Optimization A Fast Algorithm for Training Support Vector Machines》

以机器学习实战数据为例，编程实践SVM，来实现分类问题。

1.SVM待求解的问题：

以上为引入松弛变量的软间隔SVM对偶问题。

拉格朗日乘子法求偏导有：

KKT条件有：

最终线性模型有：

求出所有支持向量S 可以求得b

求解目标即为alpha参数。

2.SMO解法

由以上可知，最终求解问题转化为求解alpha向量，SMO解法每次迭代选择两个alpha进行优化处理，其余alpha值看做常量，又由于约束

因此可以转化为二元在限定区间的极值问题，可以快速求解。主要待解决的问题有：

1.如何选取两个alpha.

2.选取alpha一定，如何求解alpha的最优值

先看第二个问题，以论文符号为例，对于软间隔对偶问题形式变换如下：

其中K为核函数。

第i个数据样本经过SVM的预测值为：

根据KKT有：

假设选取的参数变量为a1,a2，则目标函数可以化简为：

其中

即优化目标为：

约束条件：

由于yi取值为-1或者+1，约束条件两边同乘以y1有： yi*yi=1

a1+a2*y1*y2= $\varepsilon$ *y1 即： s=y1*y2

代入目标优化函数，消去y1转化为关于a2的二次函数：

求导得极值点：

即有：

即a2的递推式其中y2为该样本的真实类别

E为错误差距值

又由于约束条件：因此目标函数的极值并不一定在极值点取得，a2的取值有一定的范围空间。

关于a2的参数范围的确定：

约束条件为：

如下图所示，纵轴代表a2, yi取值为-1或者+1 当y1==y2时以及y1!=y2时，有以下情况

当y1!=y2时，有a1-a2=k

当0<k<c时直线为左边右下确定的a2范围为[-k,c-k]

当-c<k<0时，直线为左边左上确定的a2范围为[-k,c] 又结合0<ai<C

有 [max(0,-k),min(c,c-k)]

当y1==y2时有，a1+a2=k

0<k<c时，直线右边左下 [0，k]

k>c时，直线右边右上 [k-c,c]

有 [max(0,k-c),min(k,c)]

根据C,a1,a2关系求出a2的变化范围取值[L,H]，结合前面二次函数求极值有：

由条件

令s=y1*y2有：

以上便是a2,a1的求解过程，可以将公式简化为：

参数b的计算：

更新完ai，若0<ai<c，即此时样本为支持向量，此时b的值很准确。此时样本点（xi,yi）满足:

消去其余部分有：

每次更新两个ai,若ai对应都是支持向量，则取均值即可。

第一个问题，a1,a2如何选取：

采用启发式算法寻找参数：

寻找和时，找那些违反KKT条件的。

第一个参数：

①遍历一遍整个数据集，对每个不满足KKT条件的参数，选作第一个待修改参数

②在上面对整个数据集遍历一遍后，选择那些参数满足的子集，开始遍历，如果发现一个不满足KKT条件的，作为第一个待修改参数，然后找到第二个待修改的参数并修改，修改完后，重新开始遍历这个子集

③遍历完子集后，重新开始①②，直到在执行①和②时没有任何修改就结束

第二个参数：

①启发式找，找能让最大的

②寻找一个随机位置的满足的可以优化的参数进行修改

③在整个数据集上寻找一个随机位置的可以优化的参数进行修改

④都不行那就找下一个第一个参数

3.数据集准备

数据格式如下：

样本的特征向量是二维类别信息为-1或者+1

def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

如上所示：加载数据，前两列作为样本特征矩阵返回，最后一列作为类别列向量返回。

3.简化版SMO算法实现

每次迭代选择两个alpha进行优化处理。其中一个增加另外一个减小。同时应满足这两个alpha在间隔边界之外，并且没有进行过区间化处理或者不在边界上。

简化版SMO应用于小量数据集：在数据集上遍历每一个alpha，若alpha可以优化, 则在剩余中随机选择一个alpha,组成一个alpha对。

伪代码如下：

创建一个alpha向量并初始化为0向量

当迭代次数小于最大迭代次数时

对数据集中的每个数据向量

如果该数据向量可以被优化

随机选择另外一个数据向量

同时优化这两个数据向量

如果两个都不能被优化，退出改层循环

如果所有向量都没被优化，继续进行下一次迭代

代码如下：

随机返回0-m之间的不等于i的随机数

#随机返回一个0-m之间并且不等于i的数字
def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

二次优化问题在一定范围的极值点，要么边界处，要么极点处

#根据边界条件以及二次优化求得极值问题，即极值要么在边界处取得，要么在极值点处取得
def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

简化版SMO算法：输入为数据集，类别集，参数C, 容错率，最大循环次数返回参数alpha以及b

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
	#数据矩阵  类别列向量
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
	#参数b  以及样本总数  特征总数
    b = 0; m,n = shape(dataMatrix)
	#alpha向量的初始化
    alphas = mat(zeros((m,1)))
	#迭代次数
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
		#遍历每一个数据集
        for i in range(m):
			#样本xi的预测值
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
			#误差
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
			#如果违反KKT 并且误差不在允许的范围内  可以进行优化
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
				#随机选择第二个参数
                j = selectJrand(i,m)
				#第二个样本的预测值
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
				#第二个样本的误差
                Ej = fXj - float(labelMat[j])
				#保存参数的旧值
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
				#求出alpha2的取值范围
                if (labelMat[i] != labelMat[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 L==H: print ("L==H"); continue
				#
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print ("eta>=0"); continue
				#第二个参数进行迭代
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
				#结合边界条件 确定极值点
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); continue
				#第一个参数进行改变
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
			    #进行b1,b2的计算
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
				#判断哪个对应支持向量  若没有则取均值
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print ("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print ("iteration number: %d" % iter)
    return b,alphas

1.预测值的计算：

将求和转换为矩阵运算有：

fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b

2.判断是否进行优化的条件：

(labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)

该alpha满足0<alpha<c,并且误差不在允许的范围内。

3.进行alpha2参数的范围判断：

if (labelMat[i] != labelMat[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])

根据两个样本的类别信息，进行参数alpha2的范围判断

y1!=y2

y1==y2

4.进行参数alpha2的调整

eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
				#第二个参数进行迭代
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
				#结合边界条件 确定极值点
                alphas[j] = clipAlpha(alphas[j],H,L)

具体解释由以下三式给出：

5.参数alpha1的计算：

alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])

其中 s=y1*y2

6.参数b的计算：

b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
				#判断哪个对应支持向量  若没有则取均值
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0

b的迭代式如上所示：观察alpha1以及alpha2是否为支持向量，若都不是则取均值作为b.

运行显示如下：

显示参数向量alpha中大于0的部分

输出支持向量：

由，可以计算w

def calcWs(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

作图绘出数据点分布以及线性平面，支持向量：

def plotfig_SVM(xMat,yMat,ws,b,alphas):
    from matplotlib import pyplot as plt
    xMat = mat(xMat)
    yMat = mat(yMat)
    b = array(b)[0] #b原来是矩阵，先转为数组类型后其数组大小为（1,1），所以后面加[0]，变为(1,)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xMat[:,0].flatten().A[0],xMat[:,1].flatten().A[0]) #注意flatten的用法
    x = arange(-1.0,10.0,0.1) #x最大值，最小值根据原数据集dataArr[:,0]的大小而定
    y =(-b-ws[0][0]*x)/ws[1][0] #根据x.w + b = 0 得到，其式子展开为w0.x1 + w1.x2 + b = 0,x2就是y值
    ax.plot(x,y)
    for i in range(100): #找到支持向量，并在图中标红
        if alphas[i]>0.0:
            ax.plot(xMat[i,0],xMat[i,1],'ro')
    plt.show()

结果如下所示：