支持向量机 SVM 算法推导优缺点代码实现 in Python

最新推荐文章于 2025-04-24 21:59:08 发布

置顶 TangowL

最新推荐文章于 2025-04-24 21:59:08 发布

阅读量4.1k

点赞数 9

分类专栏：机器学习

本文链接：https://blog.csdn.net/lipengcn/article/details/49721251

版权

机器学习专栏收录该内容

28 篇文章

订阅专栏

1、基本思想

前面讲到的Logistic Regression在拟合过程，实际上关注所有样本点的贡献，即寻找这么一个超平面，使得正例的特征远大于0，负例的特征远小于0，强调在全部训练数据上达到这一目标。而在SVM中实际也是在拟合这么一个超平面，只是关注点集中于使得距离此分界面最近的几个点据此平面的距离最大化，可见相较Logistic Regression更关注局部，但也更实际，显然这样的分割显得更科学。

SVM中定义了新的概念——函数间隔以及几何间隔——用以描述我们分类预测的正确性和确信度。

因此，SVM的基本模型是定义在特征空间上的间隔最大的线性分类器，学习策略就是间隔最大化，可形式化为求解凸二次规划的优化问题，学习算法即求解改凸二次规划的最优化算法。当T线性可分时，通过硬间隔最大化进行学习；当T近似线性可分时，通过软间隔最大化学习；当T线性不可分时，通过软间隔最大化和核技巧进行学习。

非线性映射是SVM方法的理论基础，SVM利用内积核函数代替向高维空间的非线性映射。

下面我们将推导SVM的学习算法，并给出主流求解算法SMO。

2、SVM学习算法

2.1 线性可分SVM

2.2 线性SVM

现实问题中，训练数据集往往是线性不可分的，即某些样本点无法满足函数间隔大于等于1的约束条件，因此，需对上述学习算法进行修改，即对每个样本点引进一个非负松弛变量，使之满足约束条件，同时对于引入的松弛变量，目标函数需要支付一个代价函数，然后按照线性可分SVM学习过程学习即可。

2.3 非线性SVM

处理分线性分类问题，需要将上述线性模型改进为非线性模型，这里需要应用核技巧。

核技巧应用到SVM的思想，即通过一个非线性变换将输入空间（欧式空间或其离散集合）映射到一个特征空间（希尔伯特空间），使得在输入空间的超曲面模型对应于特征空间中的超平面模型（SVM），如此一来，分类问题的学习通过在特征空间中求解线性SVM即可。这里基于这样的事实：线性不可分的样例集合，通过映射到高维特征空间，往往就变得线性可分了。

核函数的选用依赖领域知识，其选择的有效性需通过试验验证。上图介绍了常用核函数和非线性SVM非基本学习过程。假设原输入空间为n维，则多项式核函数将其映射到了维的特征空间，高斯核函数则将其映射到了无穷维的特征空间。

3、序列最小最优化算法SMO

3.1 坐标下降法

在讲解SMO之前，先简单介绍一种坐标下降法处理OP2问题。

对于下述问题时，

除了前面用到的梯度下降法、牛顿法，还可以使用坐标下降法，即

由于内循环是个关于单变量的函数，直接求导即可。如果内循环可以很快滴达到最优，则坐标下降法是个很高效的求极值方法。

3.2 SMO算法

SMO算法用来求解前面各种SVM遗留的OP2，SMO是目前最快的二次规划优化算法，特别针对线性SVM和数据稀疏时性能尤为优秀。

按照坐标下降法思路，我们可以首先固定除某个变量alpha_i外其余所有的变量，然后在alpha_i上求极值，循环往复。但由于OP2中约束条件1的存在使得此时alpha_i实际是个定值，使得本思路不可行。因此，考虑每次选定两个变量alpha_i和alpha_j，其余全固定，此时alpha_j可由alpha_i和其余所有固定变量表示，将其再会带到W，可见此时的W是关于alpha_i一个变量的函数了，可运用坐标下降法。

下面，我们将进行SMO算法推导

3.3 算法流程

4、深入说明

4.1 支持向量

在线性可分情况下，T中与分离超平面距离最近的样本点的实例称为支持向量，即使得成立的点，即支持向量一定在间隔边界上。显然从上述算法推导可以看出，在决定分离超平面时只有支持向量起作用，其他实例点并不起作用，可见增、删非支持向量样本对模型没有影响，使得SVM具有更好的鲁棒性。

此外，由于SVM 的最终决策函数只由少数的支持向量所确定，其计算的复杂性取决于支持向量的数目，而不是样本空间的维数，这在某种意义上避免了“维数灾难”。

支持向量的个数一般很少，并且支持向量的数目存在一个最优值。如果支持向量太少，就可能会得到一个很差的决策边界；而支持向量太多，就相当于每次都利用整个数据集进行分类，此时分类方法类似于kNN。

4.2 多分类SVM

可参考文献“C. W. Huset， A Comparison of Methods for Multiclass Support Vector Machines”

4.3 优点

由于SVM算法本身的优化目标是结构化风险最小，而不是经验风险最小，避免了过学习，因此，在小样本训练集上能够得到比其它算法更良好的泛化能力。

由于它是一个凸优化问题，因此局部最优解一定是全局最优解的优点。

包括4.1中提到的健壮性等优点

4.3 缺点

SVM算法对大规模训练样本难以实施

由于SVM是借助二次规划来求解支持向量，而求解二次规划将涉及N阶矩阵的计算（N为样本的个数），当N数目很大时该矩阵的存储和计算将耗费大量的机器内存和运算时间。针对以上问题的主要改进有

J.Platt的SMO算法

T.Joachims的SVM

C.J.C.Burges等的PCGC

张学工的CSVM

以及O.L.Mangasarian等的SOR算法

对非线性问题没有通用解决方案

必须谨慎选择核函数来处理，而这个过程只能不断试验

5、Python代码实现

from numpy import *

### helper functions ###
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

def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj


### simplified SMO ###
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    N, n = shape(dataMatrix)
    alphas = mat(zeros((N, 1)))
    b = 0
    iter = 0
    while(iter < maxIter):
        alphaPairsChanged = 0 # sign to denote whether the whole alphas have been optimized
        for ii in range(N): # outter alpha, alpha_1
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix*dataMatrix[ii,:].T)) + b
            Ei = fXi - float(labelMat[ii])
            # denoting the error is large which is worth to optimized
            if ((labelMat[ii]*Ei<-toler) and (alphas[ii]<C)) or ((labelMat[ii]*Ei>toler) and (alphas[ii]>0)):
                jj = selectJrand(ii,N) # inner alpha, alpha_2
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix*dataMatrix[jj,:].T)) + b
                Ej = fXi - float(labelMat[jj])
                alphaIold = alphas[ii].copy()
                alphaJold = alphas[jj].copy()
                if (labelMat[ii] != labelMat[jj]):
                    L = max(0, alphas[jj]-alphas[ii])
                    H = min(C, C+alphas[jj]-alphas[ii])
                else:
                    L = max(0, alphas[jj]+alphas[ii]-C)
                    H = min(C, alphas[jj]+alphas[ii])
                if L == H:
                    print "L == H"
                    continue
                eta = 2.0*dataMatrix[ii,:]*dataMatrix[jj,:].T - dataMatrix[ii,:]*dataMatrix[ii,:].T - dataMatrix[jj,:]*dataMatrix[jj,:].T
                if eta >= 0:
                    print "eta >= 0"
                    continue
                alphas[jj] -= labelMat[jj]*(Ei-Ej)/eta
                alphas[jj] = clipAlpha(alphas[jj], H, L)
                if (abs(alphas[jj] - alphaJold) < 0.00001):
                    print "j not moving enough"
                    continue
                alphas[ii] += labelMat[jj]*labelMat[ii]*(alphaJold-alphas[jj])
                b1 = b - Ei - labelMat[ii]*(alphas[ii]-alphaIold)*dataMatrix[ii,:]*dataMatrix[ii,:].T - labelMat[jj]*(alphas[jj]-alphaJold)*dataMatrix[ii,:]*dataMatrix[jj,:].T
                b2 = b - Ej - labelMat[ii]*(alphas[ii]-alphaIold)*dataMatrix[ii,:]*dataMatrix[jj,:].T - labelMat[jj]*(alphas[jj]-alphaJold)*dataMatrix[jj,:]*dataMatrix[jj,:].T
                if (alphas[ii] > 0) and (alphas[ii] < C):
                    b = b1
                elif (alphas[jj] > 0) and (alphas[jj] < C):
                    b = b2
                else:
                    b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print "iter: %d i: %d, pairs changed %d" % (iter, ii, alphaPairsChanged)
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
            print "ieration number: %d" % iter
    return b, alphas
    

### helper functions for full SMO ###
def kernelTrans(X, A, kTup):
	N, n = shape(X)
	K = mat(zeros((N, 1)))
	if kTup[0] == 'lin': # linear kernel function
		K = Z*A.T
	elif kTup[0] == 'rbf': # radial basis function
		for jj in range(N):
			deltaRow = X[j,:] - A
			K[jj] = deltaRow*deltaRow.T
		K = exp(K/(-2*kTup[1]**2))
	else:
		raise NameError('Houston We Have a Problem - - That Kernel is not recognized')
	return K

class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.N = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.N, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.N, 2)))
		self.K = mat(zeros((self.N, self.N)))
		for ii in range(self.N):
			self.K[:,ii] = kernelTrans(self.X, self.X[ii,:], kTup)

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

def selectJ(i, oS, Ei):
    # select the inner alpha with the alpha-j to maximize the abs(Ei-Ej)
    maxK = -1
    maxDeltaE = 0   # abs(Ei-Ej)
    Ej = 0
    oS.eCache[i] = [1, 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):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.N)
        Ej = calcEk(oS, j)
    return j, Ej

def updataEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]
    

### searching the optimization boundary for full SMO ###
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei<-oS.tol) and (oS.alphas[i]<oS.C)) or ((oS.labelMat[i]*Ei>oS.tol) and (oS.alphas[i]>0)):
        j, Ej = selectJ(i, oS, Ei)
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        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
        eta = 2.0*oS.K[i,:]*oS.K[j,:].T - oS.K[i,:]*oS.K[i,:].T - oS.K[j,:]*oS.K[j,:].T
        if eta >= 0:
            print "eta >= 0"
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updataEk(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])
        b1 = 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]
        b2 = oS.b - Ej - oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,j]
        if (oS.alphas[i] > 0) and (oS.alphas[i] < oS.C):
            oS.b = b1
        elif (oS.alphas[j] > 0) and (oS.alphas[j] < oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0


### outter iteration for full SMO ###
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter<maxIter) and ((alphaPairsChanged>0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.N):
                alphaPairsChanged += innerL(i, oS)
				print "fullSet, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged)
            iter += 1
        else:
            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
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True
        print "iteration number: %d" % iter
    return oS.b, oS.alphas


### classify ###
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    N, n = shape(X)
    w = zeros((n, 1))
    for ii in range(N):
        w += multiply(alphas[ii]*labelMat[ii], X[ii,:].T)
    return w