Python3：《机器学习实战》之支持向量机（4）核函数及其实现

Python3：《机器学习实战》之支持向量机（4）核函数及其实现

---

• 转载请注明作者和出处：http://blog.csdn.net/u011475210
• 代码地址：https://github.com/WordZzzz/ML/tree/master/Ch06
• 操作系统：WINDOWS 10
• 软件版本：python-3.6.2-amd64
• 编  者：WordZzzz

---

前言：

  前面三篇讲的都是SVM线性分类器，如果数据集非线性可分（如下图所示），那就需要做些修改了。

  显而易见，在该数据中存在某种可以识别的模式。接下来，我们就需要使用一种称为核函数（kernel）的工具将数据转换成易于分类器理解的形式。在本篇博文中，我们先对核函数进行简单的了解，然后重点研究径向基函数（radial basis function，最流行的核函数）及其实现，最后在我们之前的手写数字识别问题上进行实践。

核函数简介

  再看本文章开头的图片，我们似乎可以用一个圆来吧数据划分开来，但是对于线性分类器来说，这好像很难实现。我们或许可以对数据进行某种形式的转换，从而得到某些新的变量来表示数据。在这个例子中，我们将数据从一个特征空间转换到另一个特征空间，在新空间下，我们可以很容易利用已有的工具对数据进行处理。这个过程，学过数学的都知道，即从一个特征空间到另一个特征空间的映射。通常情况下，核函数实现的是低维到高维的映射，以后其他算法涉及到的PCA等，则是高维到低微的映射。经过空间转换之后，我们可以在高维空间解决线性问题，这也就等价于在低维空间中解决非线性问题。

  SVM优化中一个特别好的地方就是，所有的运算都可以写成内积（inner product，也叫点积）的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以把内积运算替换成核函数，而不必做简化处理。将内积替换成核函数的方式被称为核技巧（kernel trick）或者核“变电”（kernel substation）。

  当然，核函数并不仅仅应用于SVM中，很多其他的机器学习算法也都用到核函数。接下来，我们就来介绍一个流行的核函数，那就是径向基函数。

径向基函数

  径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以使从<0,0>向量或者其他向量开始计算的距离。我们用到的径向基函数的高斯版本公式为：

$$k(x,y) = exp(\frac{{-\begin{Vmatrix} x-y \end{Vmatrix}}^2}{2 \sigma^2})$$

  其中，σ是用户定义的用于确定到达率（reach）或者说是函数值跌落到0的速度参数。

  上述高斯核函数将数据从其特征空间映射到跟高维的空间，具体来说这里是映射到一个无穷维的空间。在该数据集上，使用高斯核函数得到一个很好的结果，当然，该函数也可以用于许多其他的数据集，并且也能够得到低错误率的结果。

核函数实现

  如果在svmMLiA.py文件中添加一个函数并稍作修改，那么我们就能在已有代码中使用核函数了（所有与核函数实现相关的函数，函数名末尾都是K）。其中主要区分代码在innerLK()和calcEk()中，我已经重点标记出。

  新添加的核转换函数，主要用于填充结构体和后续的计算：

def kernelTrans(X, A, kTup):
    """
    Function：   核转换函数

    Input：      X：数据集
                A：某一行数据
                kTup：核函数信息

    Output： K：计算出的核向量
    """ 
    #获取数据集行列数
    m, n = shape(X)
    #初始化列向量
    K = mat(zeros((m, 1)))
    #根据键值选择相应核函数
    #lin表示的是线性核函数
    if kTup[0] == 'lin': K = X * A.T
    #rbf表示径向基核函数
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow * deltaRow.T
        #对矩阵元素展开计算，而不像在MATLAB中一样计算矩阵的逆
        K =  exp(K/(-1*kTup[1]**2))
    #如果无法识别，就报错
    else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    #返回计算出的核向量
    return K

  其他函数：

class optStructK:
    """
    Function：   存放运算中重要的值

    Input：      dataMatIn：数据集
                classLabels：类别标签
                C：常数C
                toler：容错率
                kTup：速度参数

    Output： X：数据集
                labelMat：类别标签
                C：常数C
                tol：容错率
                m：数据集行数
                b：常数项
                alphas：alphas矩阵
                eCache：误差缓存
                K：核函数矩阵
    """ 
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))

        """ 主要区分 """
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
        """ 主要区分 """

def calcEkK(oS, k):
    """
    Function：   计算误差值E

    Input：      oS：数据结构
                k：下标

    Output： Ek：计算的E值
    """ 

    """ 主要区分 """
    #计算fXk，整个对应输出公式f(x)=w`x + b
    #fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k,:].T)) + oS.b
    fXk = float(multiply(oS.alphas, oS.labelMat).T*oS.K[:, k] + oS.b)
    """ 主要区分 """

    #计算E值
    Ek = fXk - float(oS.labelMat[k])
    #返回计算的误差值E
    return Ek

def selectJK(i, oS, Ei):
    """
    Function：   选择第二个alpha的值

    Input：      i：第一个alpha的下标
                oS：数据结构
                Ei：计算出的第一个alpha的误差值

    Output： j：第二个alpha的下标
                Ej：计算出的第二个alpha的误差值
    """ 
    #初始化参数值
    maxK = -1; maxDeltaE = 0; Ej = 0
    #构建误差缓存
    oS.eCache[i] = [1, Ei]
    #构建一个非零列表，返回值是第一个非零E所对应的alpha值，而不是E本身
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    #如果列表长度大于1，说明不是第一次循环
    if (len(validEcacheList)) > 1:
        #遍历列表中所有元素
        for k in validEcacheList:
            #如果是第一个alpha的下标，就跳出本次循环
            if k == i: continue
            #计算k下标对应的误差值
            Ek = calcEkK(oS, k)
            #取两个alpha误差值的差值的绝对值
            deltaE = abs(Ei - Ek)
            #最大值更新
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        #返回最大差值的下标maxK和误差值Ej
        return maxK, Ej
    #如果是第一次循环，则随机选择alpha，然后计算误差
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEkK(oS, j)
    #返回下标j和其对应的误差Ej
    return j, Ej

def updateEkK(oS, k):
    """
    Function：   更新误差缓存

    Input：      oS：数据结构
                j：alpha的下标

    Output： 无
    """ 
    #计算下表为k的参数的误差
    Ek = calcEkK(oS, k)
    #将误差放入缓存
    oS.eCache[k] = [1, Ek]

def innerLK(i, oS):
    """
    Function：   完整SMO算法中的优化例程

    Input：      oS：数据结构
                i：alpha的下标

    Output： 无
    """ 
    #计算误差
    Ei = calcEkK(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)):
        #启发式选择第二个alpha值
        j, Ej = selectJK(i, oS, Ei)
        #利用copy存储刚才的计算值，便于后期比较
        alphaIold = oS.alphas[i].copy(); alpahJold = oS.alphas[j].copy();
        #保证alpha在0和C之间
        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.X[i, :]*oS.X[j, :].T - oS.X[i, :]*oS.X[i, :].T - oS.X[j, :]*oS.X[j, :].T
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
        """ 主要区分 """

        #如果最优修改量大于0，则不做处理直接跳出本次循环，这里对真实SMO做了简化处理
        if eta >= 0: print("eta>=0"); return 0
        #计算新的alphas[j]的值
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        #对新的alphas[j]进行阈值处理
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        #更新误差缓存
        updateEkK(oS, j)
        #如果新旧值差很小，则不做处理跳出本次循环
        if (abs(oS.alphas[j] - alpahJold) < 0.00001): print("j not moving enough"); return 0
        #对i进行修改，修改量相同，但是方向相反
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alpahJold - oS.alphas[j])
        #更新误差缓存
        updateEkK(oS, i)

        """ 主要区分 """
        #更新常数项
        #b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :]*oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * oS.X[i, :]*oS.X[j, :].T
        #b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :]*oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * oS.X[j, :]*oS.X[j, :].T
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * 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] - alpahJold) * oS.K[j, j]
        """ 主要区分 """

        #谁在0到C之间，就听谁的，否则就取平均值
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[i]): oS.b = b2
        else: oS.b = (b1 + b2) / 2.0
        #成功返回1
        return 1
    #失败返回0
    else: return 0

def smoPK(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin', 0)):
    """
    Function：   完整SMO算法

    Input：      dataMatIn：数据集
                classLabels：类别标签
                C：常数C
                toler：容错率
                maxIter：最大的循环次数
                kTup：速度参数

    Output： b：常数项
                alphas：数据向量
    """ 
    #新建数据结构对象
    oS = optStructK(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    #初始化迭代次数
    iter = 0
    #初始化标志位
    entireSet = True; alphaPairsChanged = 0
    #终止条件：迭代次数超限、遍历整个集合都未对alpha进行修改
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        #根据标志位选择不同的遍历方式
        if entireSet:
            #遍历任意可能的alpha值
            for i in range(oS.m):
                #选择第二个alpha值，并在可能时对其进行优化处理
                alphaPairsChanged += innerLK(i, oS)
                print("fullSet, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
            #迭代次数累加
            iter += 1
        else:
            #得出所有的非边界alpha值
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            #遍历所有的非边界alpha值
            for i in nonBoundIs:
                #选择第二个alpha值，并在可能时对其进行优化处理
                alphaPairsChanged += innerLK(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

  接下来我们写测试函数。整个代码中最重要的是for循环开始的那两行，他们给出了如何利用核函数进行分类。首先利用结构初始化方法中使用过的kernelTrans()函数，得到转换后的数据。然后，再用其与前面的alpha及类别标签值求积。特别需要注意观察的是，我们是如何做到只需要支持向量数据就可以进行分类的。

def testRbf(k1 = 1.3):
    """
    Function：   利用核函数进行分类的径向基测试函数

    Input：      k1：径向基函数的速度参数

    Output： 输出打印信息
    """ 
    #导入数据集
    dataArr, labelArr = loadDataSet('testSetRBF.txt')
    #调用Platt SMO算法
    b, alphas = smoPK(dataArr, labelArr, 200, 0.00001, 10000, ('rbf', k1))
    #初始化数据矩阵和标签向量
    datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
    #记录支持向量序号
    svInd = nonzero(alphas.A > 0)[0]
    #读取支持向量
    sVs = datMat[svInd]
    #读取支持向量对应标签
    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))

  上述代码分别在训练集和测试集上进行性能测试，打印输出如下：

>>> reload(svmMLiA)
<module 'svmMLiA' from 'E:\\机器学习实战\\mycode\\Ch06\\svmMLiA.py'>
>>> svmMLiA.testRbf()
L==H
fullSet, iter: 0 i: 0, pairs changed 0
fullSet, iter: 0 i: 1, pairs changed 1
fullSet, iter: 0 i: 2, pairs changed 2
fullSet, iter: 0 i: 3, pairs changed 3
···
fullSet, iter: 6 i: 96, pairs changed 0
fullSet, iter: 6 i: 97, pairs changed 0
fullSet, iter: 6 i: 98, pairs changed 0
fullSet, iter: 6 i: 99, pairs changed 0
iteration number: 7
there are 27 Support Vectors
the training error rate is: 0.030000
the test error rate is: 0.040000

  当然，大家也可以尝试更换不同的k1参数以观察测试错误率、训练错误率、支持向量个数随k1的变化情况。下面个两张图是书上的举例。

  我们会发现，支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量太少，就可能会得到一个很差的决策边界；如果支持向量太多，也就是相当于每次都利用整个数据集进行分类，这种分类方法成为kNN（多么熟悉）。

手写识别问题回顾

  SVM对kNN的优点在于，SVM只需要保留支持向量就可以获得可比的效果，占用内存大大减小。下面，我们就用SVM来对手写数字进行分类识别（不加修改的SVM只能用于二分类问题，在这里，我们规定，如果是数字9，则为-1，否则为+1）。

def img2vector(filename):
    """
    Function：   32*32图像转换为1*1024向量

    Input：      filename：文件名称字符串

    Output： returnVect：转换之后的1*1024向量
    """ 
    #初始化要返回的1*1024向量
    returnVect = zeros((1, 1024))
    #打开文件
    fr = open(filename)
    #读取文件信息
    for i in range(32):
        #循环读取文件的前32行
        lineStr = fr.readline()
        for j in range(32):
            #将每行的头32个字符存储到要返回的向量中
            returnVect[0, 32*i+j] = int(lineStr[j])
    #返回要输出的1*1024向量
    return returnVect

def loadImages(dirName):
    """
    Function：   加载图片

    Input：      dirName：文件路径

    Output： trainingMat：训练数据集
                hwLabels：数据标签
    """ 
    from os import listdir
    #初始化数据标签
    hwLabels = []
    #读取文件列表
    trainingFileList = listdir(dirName)
    #读取文件个数
    m = len(trainingFileList)
    #初始化训练数据集
    trainingMat = zeros((m,1024))
    #填充数据集
    for i in range(m):
        #遍历所有文件
        fileNameStr = trainingFileList[i]
        #提取文件名称
        fileStr = fileNameStr.split('.')[0]
        #提取数字标识
        classNumStr = int(fileStr.split('_')[0])
        #数字9记为-1类
        if classNumStr == 9: hwLabels.append(-1)
        #其他数字记为+1类
        else: hwLabels.append(1)
        #提取图像向量，填充数据集
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
    #返回数据集和数据标签
    return trainingMat, hwLabels

def testDigits(kTup = ('rbf',10)):
    """
    Function：   手写数字分类函数

    Input：      kTup：核函数采用径向基函数

    Output： 输出打印信息
    """ 
    #导入数据集
    dataArr, labelArr = loadImages('trainingDigits')
    #调用Platt SMO算法
    b, alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    #初始化数据矩阵和标签向量
    datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
    #记录支持向量序号
    svInd = nonzero(alphas.A > 0)[0]
    #读取支持向量
    sVs = datMat[svInd]
    #读取支持向量对应标签
    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,:],kTup)
        #仅用支持向量预测分类
        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 = loadImages('testDigits')
    #初始化误差计数
    errorCount = 0
    #初始化数据矩阵和标签向量
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    #获取数据集行列值
    m,n = shape(datMat)
    #遍历每一行，利用核函数对测试集进行分类
    for i in range(m):
        #利用核函数转换数据
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        #仅用支持向量预测分类
        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) 

  上面的大部分代码我们都已经很熟悉了，这里就不再赘述，直接进行测试。

>>> reload(svmMLiA)
>>> svmMLiA.testDigits(('rbf', 20))
L==H
fullSet, iter: 6 i: 398, pairs changed 0
j not moving enough
fullSet, iter: 6 i: 399, pairs changed 0
fullSet, iter: 6 i: 400, pairs changed 0
j not moving enough
fullSet, iter: 6 i: 401, pairs changed 0
iteration number: 7
there are 58 Support Vectors
the training error rate is: 0.000000
the test error rate is: 0.021505

  尝试不同的σ，并尝试了线性核函数，总结得到如图结果：

Kernel,settings	Training error (%)	Test error (%)	Support vectors
RBF,0.1	0	52	402
RBF,5	0	3.2	402
RBF,10	0	0.5	99
RBF,50	0.2	2.2	41
RBF,100	1.5	4.3	26
Linear	2.7	2.2	28

  结果表明，当径向基函数的参数σ取10左右时，就可以得到最小的测试错误率。同时，最小的训练错误率，并不对应于最小的支持向量数目。另外，线性核函数的效果并不是特别糟糕，可以以牺牲线性核函数的错误率来换取分类速度的提高。

总结

  支持向量机是一种分类器。之所以称为“机”是因为它会产生一个二值决策结果，即它是一种决策“机”。支持向量机的泛化错误率较低，具有良好的学习能力，且学到的结果具有很好的推广性。这些优点使得支持向量机十分流行，有些人认为他是监督学习中最好的定式算法。

系列教程持续发布中，欢迎订阅、关注、收藏、评论、点赞哦～～(￣▽￣～)～

完的汪(∪｡∪)｡｡｡zzz